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
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
|
<?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 - 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="quad.xhtml"> Quad-precision 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="#faq">Frequently asked questions</a></li>
<li><a href="#vectorizing">Vectorizing calls to scalar functions</a></li>
<!-- <li><a href="#gnuabi">About the GNUABI version of the library</a></li>-->
<li><a href="#lto">Using link time optimization</a></li>
<li><a href="#inline">Using header files of inlinable functions</a></li>
<li><a href="#wasm">Utilizing SLEEF for WebAssembly</a></li>
<li><a href="#dispatcher">How the dispatcher works</a></li>
<li><a href="#ulp">ULP, gradual underflow and flush-to-zero mode</a></li>
<li><a href="#paynehanek">Explanatory source code for the modified Payne Hanek reduction method</a></li>
<li><a href="#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 computes directly with 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 cases 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 <a class="underlined"
href="#paynehanek">vectorized version of Payne Hanek range
reduction</a>, and all the trigonometric functions return a correct
value with the specified accuracy.
</p>
<br/>
<br/>
<p class="noindent">
<b>Q3:</b> What can I do to make sleef run faster?
</p>
<br/>
<p class="noindent">
<b>A3:</b> The most important thing is to choose the fastest
available vector extension. SLEEF is optimized for computers with
<a class="underlined"
href="https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation#Fused_multiply%E2%80%93add">FMA</a>
instructions, and it runs slow on Ivy Bridge or older CPUs and Atom,
that do not have FMA instructions. If you are not sure, use the
dispatcher. <a class="underlined" href="#dispatcher">The dispatcher in
SLEEF is not slow</a>. If you want to further speed up computation,
try using <a class="underlined" href="#lto">LTO</a>. By using LTO, the
compiler fuses the code within the library to the code calling the
library functions, and this sometimes results in considerable
performance boost. In this case, you should not use the dispatcher,
and you should use the same compiler with the same version to build
SLEEF and the program against which SLEEF is linked.
</p>
<h2 id="vectorizing">Vectorizing calls to scalar functions</h2>
<p class="noindent">
Recent x86_64 gcc can <a class="underlined"
href="https://gcc.gnu.org/projects/tree-ssa/vectorization.html">auto-vectorize</a>
calls to functions. In order to utilize this functionality, OpenMP
SIMD pragmas can be added to declarations of scalar functions
like <b class="func">Sleef_sin_u10</b> by defining
<b>SLEEF_ENABLE_OMP_SIMD</b> macro before including <b>sleef.h</b> on
x86_64 computers. With these pragmas, gcc can use its auto-vectorizer
to vectorize calls to these scalar functions. For
example, <a class="underlined" href="sophomore.c">the following
code</a> can be vectorized by gcc-10.
</p>
<pre class="code">
<code>#include <stdio.h></code>
<code></code>
<code>#define SLEEF_ENABLE_OMP_SIMD</code>
<code>#include "sleef.h"</code>
<code></code>
<code>#define N 65536</code>
<code>#define M (N + 3)</code>
<code></code>
<code>static double func(double x) { return Sleef_pow_u10(x, -x); }</code>
<code></code>
<code>double int_simpson(double a, double b) {</code>
<code> double h = (b - a) / M;</code>
<code> double sum_odd = 0.0, sum_even = 0.0;</code>
<code> for(int i = 1;i <= M-3;i += 2) {</code>
<code> sum_odd += func(a + h * i);</code>
<code> sum_even += func(a + h * (i + 1));</code>
<code> }</code>
<code> return h / 3 * (func(a) + 4 * sum_odd + 2 * sum_even + 4 * func(b - h) + func(b));</code>
<code>}</code>
<code></code>
<code>int main() {</code>
<code> double sum = 0;</code>
<code> for(int i=1;i<N;i++) sum += Sleef_pow_u10(i, -i);</code>
<code> printf("%g %g\n", int_simpson(0, 1), sum);</code>
<code>}</code>
</pre>
<pre class="command">$ gcc-10 -fopenmp -ffast-math -mavx2 -O3 sophomore.c -lsleef -S -o- | grep _ZGV
call _ZGVdN4vv_Sleef_pow_u10@PLT
call _ZGVdN4vv_Sleef_pow_u10@PLT
call _ZGVdN4vv_Sleef_pow_u10@PLT
call _ZGVdN4vv_Sleef_pow_u10@PLT
call _ZGVdN4vv_Sleef_pow_u10@PLT
call _ZGVdN4vv_Sleef_pow_u10@PLT
call _ZGVdN4vv_Sleef_pow_u10@PLT
$ █
</pre>
<!--
<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>. The auto-vectorizer in x86_64 gcc is capable of
vectorizing calls to the standard math functions and generates calls
to <a class="underlined"
href="https://sourceware.org/glibc/wiki/libmvec">libmvec</a>. The
GNUABI version of SLEEF library can be used as a substitute for
libmvec.
</p>
<pre class="code">
<code>#include <stdio.h></code>
<code>#include <math.h></code>
<code></code>
<code>#define N 65536</code>
<code>#define M (N + 3)</code>
<code></code>
<code>static double func(double x) { return pow(x, -x); }</code>
<code></code>
<code>double int_simpson(double a, double b) {</code>
<code> double h = (b - a) / M;</code>
<code> double sum_odd = 0.0, sum_even = 0.0;</code>
<code> for(int i = 1;i <= M-3;i += 2) {</code>
<code> sum_odd += func(a + h * i);</code>
<code> sum_even += func(a + h * (i + 1));</code>
<code> }</code>
<code> return h / 3 * (func(a) + 4 * sum_odd + 2 * sum_even + 4 * func(b - h) + func(b));</code>
<code>}</code>
<code></code>
<code>int main() {</code>
<code> double sum = 0;</code>
<code> for(int i=1;i<N;i++) sum += pow(i, -i);</code>
<code> printf("%g %g\n", int_simpson(0, 1), sum);</code>
<code>}</code>
</pre>
<p>
For example, <a class="underlined" href="sophomore2.c">the above
code</a> can be linked against libsleefgnuabi as shown below. You
have to specify <b>-lsleefgnuabi</b> compiler option
before <b>-lm</b> option.
</p>
<pre class="command">$ gcc-10 -ffast-math -O3 sophomore2.c -lsleefgnuabi -lm -L./lib
$ ldd a.out
linux-vdso.so.1 (0x00007ffd0c5ff000)
libsleefgnuabi.so.3 => not found
libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007f73f5f98000)
libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f73f5da6000)
/lib64/ld-linux-x86-64.so.2 (0x00007f73f60fb000)
$ LD_LIBRARY_PATH=./lib ./a.out
1.29127 1.29129
$ █
</pre>
-->
<h2 id="lto">Using link time optimization</h2>
<p class="noindent">
<a class="underlined"
href="https://en.wikipedia.org/wiki/Interprocedural_optimization">Link
time optimization (LTO)</a> is a functionality implemented in gcc,
clang and other compilers for optimizing across translation units (or
source files.) This can sometimes dramatically improve the
performance of the code, because it is capable of fusing library
functions into the code calling those functions. The build system in
SLEEF supports LTO and thus it can be built with LTO support by just
specifying <b>-DSLEEF_ENABLE_LTO=TRUE</b> cmake option. However, there are a
few things to note in order to get the optimal performance. 1. You
should not use the dispatcher with LTO. Dispatchers prevent the
functions from being fused with LTO. 2. You have to use the same
compiler with the same version to build the library and your
code. 3. You cannot build shared libraries with LTO.
</p>
<h2 id="inline">Using header files of inlinable functions</h2>
<p class="noindent">
Although LTO is considered to be a smart technique for improving the
performance of the library functions, there are difficulties in using
this functionality in real situations. One of the reasons is that
people still need to use old compilers to build their projects. SLEEF
can generate header files in which the library functions are all
defined as inline functions. This can be compiled with old compilers.
In theory, inline functions should give similar performance to LTO,
but in reality, inline functions are better. In order to generate
those header files, specify <b>-DSLEEF_BUILD_INLINE_HEADERS=TRUE</b> cmake
option. Below is an example code utilizing the generated header files
for SSE2 and AVX2. You cannot use a dispatcher with these header
files.
</p>
<pre class="code">
<code>#include <stdio.h></code>
<code>#include <stdint.h></code>
<code>#include <string.h></code>
<code>#include <x86intrin.h></code>
<code></code>
<code>#include <sleefinline_sse2.h></code>
<code>#include <sleefinline_avx2128.h></code>
<code></code>
<code>int main(int argc, char **argv) {</code>
<code> __m128d va = _mm_set_pd(2, 10);</code>
<code> __m128d vb = _mm_set_pd(3, 20);</code>
<code></code>
<code> __m128d vc = Sleef_powd2_u10sse2(va, vb);</code>
<code></code>
<code> double c[2];</code>
<code> _mm_storeu_pd(c, vc);</code>
<code></code>
<code> printf("%g, %gn", c[0], c[1]);</code>
<code></code>
<code> __m128d vd = Sleef_powd2_u10avx2128(va, vb);</code>
<code></code>
<code> double d[2];</code>
<code> _mm_storeu_pd(d, vd);</code>
<code></code>
<code> printf("%g, %gn", d[0], d[1]);</code>
<code>}</code>
</pre>
<pre class="command">$ gcc-10 -ffp-contract=off -O3 -march=native helloinline.c -I./include
$ ./a.out
1e+20, 8
1e+20, 8
$ nm -g a.out
00000000000036a0 R Sleef_rempitabdp
0000000000003020 R Sleef_rempitabsp
0000000000003000 R _IO_stdin_used
w _ITM_deregisterTMCloneTable
w _ITM_registerTMCloneTable
000000000000d010 D __TMC_END__
000000000000d010 B __bss_start
w __cxa_finalize@@GLIBC_2.2.5
000000000000d000 D __data_start
000000000000d008 D __dso_handle
w __gmon_start__
00000000000020a0 T __libc_csu_fini
0000000000002030 T __libc_csu_init
U __libc_start_main@@GLIBC_2.2.5
U __printf_chk@@GLIBC_2.3.4
000000000000d010 D _edata
000000000000d018 B _end
00000000000020a8 T _fini
0000000000001f40 T _start
000000000000d000 W data_start
0000000000001060 T main
$ █
</pre>
<h2 id="wasm">Utilizing SLEEF for WebAssembly</h2>
<p class="noindent">
Since <a class="underlined"
href="https://emscripten.org/">Emscripten</a> supports SSE2
intrinsics, the SSE2 inlinable function header can be used for
<a class="underlined" href="https://webassembly.org/">WebAssembly</a>.
</p>
<pre class="code">
<code>#include <stdio.h></code>
<code>#include <emmintrin.h></code>
<code></code>
<code>#include "sleefinline_sse2.h"</code>
<code></code>
<code>int main(int argc, char **argv) {</code>
<code> double a[] = {2, 10};</code>
<code> double b[] = {3, 20};</code>
<code></code>
<code> __m128d va, vb, vc;</code>
<code></code>
<code> va = _mm_loadu_pd(a);</code>
<code> vb = _mm_loadu_pd(b);</code>
<code></code>
<code> vc = Sleef_powd2_u10sse2(va, vb);</code>
<code></code>
<code> double c[2];</code>
<code></code>
<code> _mm_storeu_pd(c, vc);</code>
<code></code>
<code> printf("pow(%g, %g) = %gn", a[0], b[0], c[0]);</code>
<code> printf("pow(%g, %g) = %gn", a[1], b[1], c[1]);</code>
<code>}</code>
</pre>
<pre class="command">$ emcc -O3 -msimd128 -msse2 hellowasm.c
$ ../node-v15.7.0-linux-x64/bin/node --experimental-wasm-simd ./a.out.js
pow(2, 3) = 8
pow(10, 20) = 1e+20
$ █
</pre>
<h2 id="dispatcher">How the dispatchers work</h2>
<p class="noindent">
SLEEF implements versions of functions that are implemented with
each vector extension of the architecture. A dispatcher is a
function that dynamically selects the fastest implementatation for
the computer it runs. The dispatchers in SLEEF are designed to have
very low overhead.
</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 © 2010-2025 SLEEF Project, Naoki Shibata and contributors.<br/>
SLEEF is open-source software and is distributed under the Boost Software License, Version 1.0.
</p>
<script type="text/javascript">
var sc_project=13098265;
var sc_invisible=1;
var sc_security="518de45e";
</script>
<script type="text/javascript"
src="https://www.statcounter.com/counter/counter.js"
async="async"></script>
<noscript><div class="statcounter"><a title="Web Analytics"
href="https://statcounter.com/" target="_blank"><img
class="statcounter"
src="https://c.statcounter.com/13098265/0/518de45e/1/"
alt="Web Analytics"
referrerPolicy="no-referrer-when-downgrade"></img></a></div></noscript>
</body>
</html>
|
