1 // Explanatory source code for the modified Payne Hanek reduction
2 // http://dx.doi.org/10.1109/TPDS.2019.2960333
3
4 #include <stdio.h>
5 #include <stdlib.h>
6 #include <math.h>
7 #include <mpfr.h>
8
9 typedef struct { double x, y; } double2;
dd(double d)10 double2 dd(double d) { double2 r = { d, 0 }; return r; }
d2i(double d)11 int64_t d2i(double d) { union { double f; int64_t i; } tmp = {.f = d }; return tmp.i; }
i2d(int64_t i)12 double i2d(int64_t i) { union { double f; int64_t i; } tmp = {.i = i }; return tmp.f; }
upper(double d)13 double upper(double d) { return i2d(d2i(d) & 0xfffffffff8000000LL); }
clearlsb(double d)14 double clearlsb(double d) { return i2d(d2i(d) & 0xfffffffffffffffeLL); }
15
ddrenormalize(double2 t)16 double2 ddrenormalize(double2 t) {
17 double2 s = dd(t.x + t.y);
18 s.y = t.x - s.x + t.y;
19 return s;
20 }
21
ddadd(double2 x,double2 y)22 double2 ddadd(double2 x, double2 y) {
23 double2 r = dd(x.x + y.x);
24 double v = r.x - x.x;
25 r.y = (x.x - (r.x - v)) + (y.x - v) + (x.y + y.y);
26 return r;
27 }
28
ddmul(double x,double y)29 double2 ddmul(double x, double y) {
30 double2 r = dd(x * y);
31 r.y = fma(x, y, -r.x);
32 return r;
33 }
34
ddmul2(double2 x,double2 y)35 double2 ddmul2(double2 x, double2 y) {
36 double2 r = ddmul(x.x, y.x);
37 r.y += x.x * y.y + x.y * y.x;
38 return r;
39 }
40
41 // This function computes remainder(a, PI/2)
modifiedPayneHanek(double a)42 double2 modifiedPayneHanek(double a) {
43 double table[4];
44 int scale = fabs(a) > 1e+200 ? -128 : 0;
45 a = ldexp(a, scale);
46
47 // Table genration
48
49 mpfr_set_default_prec(2048);
50 mpfr_t pi, m;
51 mpfr_inits(pi, m, NULL);
52 mpfr_const_pi(pi, GMP_RNDN);
53
54 mpfr_d_div(m, 2, pi, GMP_RNDN);
55 mpfr_set_exp(m, mpfr_get_exp(m) + (ilogb(a) - 53 - scale));
56 mpfr_frac(m, m, GMP_RNDN);
57 mpfr_set_exp(m, mpfr_get_exp(m) - (ilogb(a) - 53));
58
59 for(int i=0;i<4;i++) {
60 table[i] = clearlsb(mpfr_get_d(m, GMP_RNDN));
61 mpfr_sub_d(m, m, table[i], GMP_RNDN);
62 }
63
64 mpfr_clears(pi, m, NULL);
65
66 // Main computation
67
68 double2 x = dd(0);
69 for(int i=0;i<4;i++) {
70 x = ddadd(x, ddmul(a, table[i]));
71 x.x = x.x - round(x.x);
72 x = ddrenormalize(x);
73 }
74
75 double2 pio2 = { 3.141592653589793*0.5, 1.2246467991473532e-16*0.5 };
76 x = ddmul2(x, pio2);
77 return fabs(a) < 0.785398163397448279 ? dd(a) : x;
78 }
79
main(int argc,char ** argv)80 int main(int argc, char **argv) {
81 double a = ldexp(6381956970095103.0, 797);
82 if (argc > 1) a = atof(argv[1]);
83 printf("a = %.20g\n", a);
84
85 //
86
87 mpfr_set_default_prec(2048);
88 mpfr_t pi, pio2, x, y, r;
89 mpfr_inits(pi, pio2, x, y, r, NULL);
90
91 mpfr_const_pi(pi, GMP_RNDN);
92 mpfr_mul_d(pio2, pi, 0.5, GMP_RNDN);
93
94 //
95
96 mpfr_set_d(x, a, GMP_RNDN);
97 mpfr_remainder(r, x, pio2, GMP_RNDN);
98
99 mpfr_printf("mpfr = %.64RNf\n", r);
100
101 //
102
103 double2 dd = modifiedPayneHanek(a);
104 mpfr_set_d(x, dd.x, GMP_RNDN);
105 mpfr_add_d(x, x, dd.y, GMP_RNDN);
106
107 mpfr_printf("dd = %.64RNf\n", x);
108
109 mpfr_sub(x, x, r, GMP_RNDN);
110 mpfr_abs(x, x, GMP_RNDN);
111 mpfr_div(x, x, r, GMP_RNDN);
112
113 double err = mpfr_get_d(x, GMP_RNDN);
114 printf("error = %g\n", err);
115 }
116