1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
|
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1 plus MathML 2.0//EN" "http://www.w3.org/Math/DTD/mathml2/xhtml-math11-f.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en-US">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
<meta name="viewport" content="width=device-width, initial-scale=0.4"/>
<meta name="google" content="notranslate" />
<link rel="canonical" href="https://sleef.org/additional.xhtml" />
<link rel="icon" href="favicon.png" />
<link rel="stylesheet" type="text/css" href="texlike.css"/>
<link rel="stylesheet" type="text/css" href="//fonts.googleapis.com/css?family=Ubuntu" />
<link rel="stylesheet" type="text/css" href="sleef.css"/>
<title>SLEEF - Additional Notes</title>
</head>
<body translate="no" class="notranslate">
<h1>SLEEF Documentation - Additional Notes</h1>
<h2>Table of contents</h2>
<ul class="none" style="font-family: arial, sansserif; padding-left: 0.5cm;">
<li><a class="underlined" href="index.xhtml">Introduction</a></li>
<li><a class="underlined" href="compile.xhtml">Compiling and installing the library</a></li>
<li><a class="underlined" href="purec.xhtml">Math library reference</a></li>
<li><a class="underlined" href="dft.xhtml">DFT library reference</a></li>
<li><a class="underlined" href="misc.xhtml">Other tools included in the package</a></li>
<li><a class="underlined" href="benchmark.xhtml">Benchmark results</a></li>
<li> </li>
<li><a class="underlined" href="additional.xhtml">Additional notes</a></li>
<ul class="disc">
<li><a href="additional.xhtml#faq">Frequently asked questions</a></li>
<li><a href="additional.xhtml#gnuabi">About the GNUABI version of the library</a></li>
<li><a href="additional.xhtml#dispatcher">How the dispatcher works</a></li>
<li><a href="additional.xhtml#ulp">ULP, gradual underflow and flush-to-zero mode</a></li>
<li><a href="additional.xhtml#paynehanek">Explanatory source code for the modified Payne Hanek reduction method</a></li>
<li><a href="additional.xhtml#logo">About the logo</a></li>
</ul>
</ul>
<h2 id="faq">Frequently asked questions</h2>
<p class="noindent">
<b>Q1:</b> Is the scalar functions in SLEEF faster than the
corresponding functions in the standard C library?
</p>
<br/>
<p class="noindent">
<b>A1:</b> No. Todays standard C libraries are very well optimized,
and there is small room for further optimization. The reason why
SLEEF is fast is that it carries out computation directly on SIMD
registers and ALUs. This is not simple as it sounds, because
conditional branches have to be eliminated in order to take full
advantage of SIMD computation. If the algorithm requires conditional
branches according to the argument, it must prepare for the case
where the elements in the input vector contain both values that
would make a branch happen and not happen. This would spoil the
advantage of SIMD computation, because each element in a vector
would require a different code path.
</p>
<br/>
<br/>
<p class="noindent">
<b>Q2:</b> Do the trigonometric functions (e.g. sin) in SLEEF return
correct values for the whole range of inputs?
</p>
<br/>
<p class="noindent">
<b>A2:</b> Yes. SLEEF does implement a vectorized version of Payne
Hanek range reduction, and all the trigonometric functions return
a correct value with the specified accuracy.
</p>
<h2 id="gnuabi">About the GNUABI version of the library</h2>
<p class="noindent">
The GNUABI version of the library (libsleefgnuabi.so) is built for
x86 and aarch64 architectectures. This library provides an API
compatible with <a class="underlined"
href="https://sourceware.org/glibc/wiki/libmvec">libmvec</a> in
glibc, and the API comforms to the <a class="underlined"
href="https://sourceware.org/glibc/wiki/libmvec?action=AttachFile&do=view&target=VectorABI.txt">x86
vector ABI</a>, <a class="underlined"
href="https://developer.arm.com/docs/101129/latest">AArch64 vector
ABI</a> and <a class="underlined"
href="https://github.com/power8-abi-doc/vector-function-abi/">Power
Vector ABI</a>. This library is built and installed by default, and
certain compilers call the functions in this library.
</p>
<h2 id="dispatcher">How the dispatchers work</h2>
<p class="noindent">
The dispatchers in SLEEF are designed to have very low overhead. This
overhead is so small and cannot be observed by microbenchmarking.
</p>
<p>
Fig. 7.1 shows a simplified code of our dispatcher. There is only
one exported function <b class="func">mainFunc</b>. When
<b class="func">mainFunc</b> is called for the first
time, <b class="func">dispatcherMain</b> is called internally,
since <i class="var">funcPtr</i> is initialized to the pointer to
<b class="func">dispatcherMain</b> (line 14). It then detects if the
CPU supports SSE 4.1 (line 7), and
rewrites <i class="var">funcPtr</i> to a pointer to the function
that utilizes SSE 4.1 or SSE 2, depending on the result of CPU
feature detection (line 10). When
<b class="func">mainFunc</b> is called for the second time, it does
not execute the
<b class="func">dispatcherMain</b>. It just executes the function
pointed by the pointer stored in <i class="var">funcPtr</i> during
the execution of
<b class="func">dispatcherMain</b>.
</p>
<p>
There are advantages in our dispatcher. The first advantage is that
it does not require any compiler-specific extension. The second
advantage is simplicity. There are only 18 lines of simple
code. Since the dispatchers are completely separated for each
function, there is not much room for bugs to get in.
</p>
<p>
The third advantage is low overhead. You might think that the
overhead is one function call including execution of the prologue
and the epilogue. However, modern compilers are smart enough to
eliminate redundant execution of the prologue, epilogue and return
instruction. The actual overhead is just one jmp instruction, which
has very small overhead since it is not conditional. This overhead
is likely hidden by out-of-order execution.
</p>
<p>
The fourth advantage is thread safety. There is only one variable
shared among threads, which is <i class="var">funcPtr</i>. There are
only two possible values for this pointer variable. The first value
is the pointer to the <b class="func">dispatcherMain</b>, and the
second value is the pointer to either <b class="func">funcSSE2</b>
or <b class="func">funcSSE4</b>, depending on the availability of
extensions. Once <i class="var">funcPtr</i> is substituted with the
pointer to <b class="func">funcSSE2</b>
or <b class="func">funcSSE4</b>, it will not be changed in the
future. It should be easy to confirm that the code works in all the
cases.
</p>
<pre class="code">
<code>static double (*funcPtr)(double arg);</code>
<code></code>
<code>static double dispatcherMain(double arg) {</code>
<code> double (*p)(double arg) = funcSSE2;</code>
<code></code>
<code>#if the compiler supports SSE4.1</code>
<code> if (SSE4.1 is available on the CPU) p = funcSSE4;</code>
<code>#endif</code>
<code></code>
<code> funcPtr = p;</code>
<code> return (*funcPtr)(arg);</code>
<code>}</code>
<code></code>
<code>static double (*funcPtr)(double arg) = dispatcherMain;</code>
<code></code>
<code>double mainFunc(double arg) {</code>
<code> return (*funcPtr)(arg);</code>
<code>}</code>
</pre>
<p style="text-align:center; margin-bottom: 1.0cm;">
Fig. 7.1: Simplified code of our dispatcher
</p>
<h2 id="ulp">ULP, gradual underflow and flush-to-zero mode</h2>
<p class="noindent">
ULP stands for "unit in the last place", which is sometimes used for
representing accuracy of calculation. 1 ULP is the distance between
the two closest floating point number, which depends on the exponent
of the FP number. The accuracy of calculation by reputable math
libraries is usually between 0.5 and 1 ULP. Here, the accuracy means
the largest error of calculation. SLEEF math library provides
multiple accuracy choices for most of the math functions. Many
functions have 3.5-ULP and 1-ULP versions, and 3.5-ULP versions are
faster than 1-ULP versions. If you care more about execution speed
than accuracy, it is advised to use the 3.5-ULP versions along with
-ffast-math or "unsafe math optimization" options for the compiler.
</p>
<p>
Note that 3.5 ULPs of error is small enough in many applications. If
you do not manage the error of computation by carefully ordering
floating point operations in your code, you would easily have that
amount of error in the computation results.
</p>
<p>
In IEEE 754 standard, underflow does not happen abruptly when the
exponent becomes zero. Instead, when a number to be represented is
smaller than a certain value, a denormal number is produced which
has less precision. This is sometimes called gradual underflow. On
some processor implementation, a flush-to-zero mode is used since it
is easier to implement by hardware. In flush-to-zero mode, numbers
smaller than the smallest normalized number are replaced with
zero. FP operations are not IEEE-754 conformant if a flush-to-zero
mode is used. A flush-to-zero mode influences the accuracy of
calculation in some cases. The smallest normalized precision number
can be referred with DBL_MIN for double precision, and FLT_MIN for
single precision. The naming of these macros is a little bit
confusing because DBL_MIN is not the smallest double precision
number.
</p>
<p>
You can see known maximum errors in math functions in glibc
at <a class="underlined"
href="http://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html">
this page.</a>
</p>
<h2 id="paynehanek">Explanatory source code for our modified Payne Hanek reduction method</h2>
<p class="noindent">
In order to evaluate a trigonometric function with a large argument,
an argument reduction method is used to find an FP remainder of
dividing the argument <i class="var">x</i> by π. We devised a
variation of the Payne-Hanek argument reduction method which is
suitable for vector computation. Fig. 7.2
shows <a class="underlined" href="ph.c">an explanatory source
code</a> for this method. See <a class="underlined"
href="http://dx.doi.org/10.1109/TPDS.2019.2960333">our paper</a> for
the details.
</p>
<pre class="code">
<code>#include <stdio.h></code>
<code>#include <stdlib.h></code>
<code>#include <math.h></code>
<code>#include <mpfr.h></code>
<code></code>
<code>typedef struct { double x, y; } double2;</code>
<code>double2 dd(double d) { double2 r = { d, 0 }; return r; }</code>
<code>int64_t d2i(double d) { union { double f; int64_t i; } tmp = {.f = d }; return tmp.i; }</code>
<code>double i2d(int64_t i) { union { double f; int64_t i; } tmp = {.i = i }; return tmp.f; }</code>
<code>double upper(double d) { return i2d(d2i(d) & 0xfffffffff8000000LL); }</code>
<code>double clearlsb(double d) { return i2d(d2i(d) & 0xfffffffffffffffeLL); }</code>
<code></code>
<code>double2 ddrenormalize(double2 t) {</code>
<code> double2 s = dd(t.x + t.y);</code>
<code> s.y = t.x - s.x + t.y;</code>
<code> return s;</code>
<code>}</code>
<code></code>
<code>double2 ddadd(double2 x, double2 y) {</code>
<code> double2 r = dd(x.x + y.x);</code>
<code> double v = r.x - x.x;</code>
<code> r.y = (x.x - (r.x - v)) + (y.x - v) + (x.y + y.y);</code>
<code> return r;</code>
<code>}</code>
<code></code>
<code>double2 ddmul(double x, double y) {</code>
<code> double2 r = dd(x * y);</code>
<code> r.y = fma(x, y, -r.x);</code>
<code> return r;</code>
<code>}</code>
<code></code>
<code>double2 ddmul2(double2 x, double2 y) {</code>
<code> double2 r = ddmul(x.x, y.x);</code>
<code> r.y += x.x * y.y + x.y * y.x;</code>
<code> return r;</code>
<code>}</code>
<code></code>
<code>// This function computes remainder(a, PI/2)</code>
<code>double2 modifiedPayneHanek(double a) {</code>
<code> double table[4];</code>
<code> int scale = fabs(a) > 1e+200 ? -128 : 0;</code>
<code> a = ldexp(a, scale);</code>
<code></code>
<code> // Table genration</code>
<code></code>
<code> mpfr_set_default_prec(2048);</code>
<code> mpfr_t pi, m;</code>
<code> mpfr_inits(pi, m, NULL);</code>
<code> mpfr_const_pi(pi, GMP_RNDN);</code>
<code></code>
<code> mpfr_d_div(m, 2, pi, GMP_RNDN);</code>
<code> mpfr_set_exp(m, mpfr_get_exp(m) + (ilogb(a) - 53 - scale));</code>
<code> mpfr_frac(m, m, GMP_RNDN);</code>
<code> mpfr_set_exp(m, mpfr_get_exp(m) - (ilogb(a) - 53));</code>
<code></code>
<code> for(int i=0;i<4;i++) {</code>
<code> table[i] = clearlsb(mpfr_get_d(m, GMP_RNDN));</code>
<code> mpfr_sub_d(m, m, table[i], GMP_RNDN);</code>
<code> }</code>
<code></code>
<code> mpfr_clears(pi, m, NULL);</code>
<code></code>
<code> // Main computation</code>
<code></code>
<code> double2 x = dd(0);</code>
<code> for(int i=0;i<4;i++) {</code>
<code> x = ddadd(x, ddmul(a, table[i]));</code>
<code> x.x = x.x - round(x.x);</code>
<code> x = ddrenormalize(x);</code>
<code> }</code>
<code></code>
<code> double2 pio2 = { 3.141592653589793*0.5, 1.2246467991473532e-16*0.5 };</code>
<code> x = ddmul2(x, pio2);</code>
<code> return fabs(a) < 0.785398163397448279 ? dd(a) : x;</code>
<code>}</code>
</pre>
<p style="text-align:center; margin-bottom: 1.0cm;">
<a href="ph.c">Fig. 7.2: Explanatory source code for our modified Payne Hanek reduction method</a>
</p>
<h2 id="logo">About the logo</h2>
<p>
It is a soup ladle. A sleef means a soup ladle in Dutch.
</p>
<br/>
<p style="text-align:center; margin-top:1cm;">
<a class="nothing" href="sleeflogo3.svg">
<img src="sleeflogo3.png" alt="logo" width="40%" height="40%" />
</a>
<br />
Fig. 7.2: SLEEF logo
</p>
<p class="footer">
Copyright © <!--YEAR--> SLEEF Project.<br/>
SLEEF is open-source software and is distributed under the Boost Software License, Version 1.0.
</p>
<!--TEST-->
</body>
</html>
|
