1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | #include "libm.h" |
21 | |
22 | #if FLT_EVAL_METHOD0==0 || FLT_EVAL_METHOD0==1 |
23 | #define EPS2.22044604925031308085e-16 DBL_EPSILON2.22044604925031308085e-16 |
24 | #elif FLT_EVAL_METHOD0==2 |
25 | #define EPS2.22044604925031308085e-16 LDBL_EPSILON1.0842021724855044340e-19L |
26 | #endif |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | static const double |
38 | toint = 1.5/EPS2.22044604925031308085e-16, |
39 | invpio2 = 6.36619772367581382433e-01, |
40 | pio2_1 = 1.57079632673412561417e+00, |
41 | pio2_1t = 6.07710050650619224932e-11, |
42 | pio2_2 = 6.07710050630396597660e-11, |
43 | pio2_2t = 2.02226624879595063154e-21, |
44 | pio2_3 = 2.02226624871116645580e-21, |
45 | pio2_3t = 8.47842766036889956997e-32; |
46 | |
47 | |
48 | int __rem_pio2(double x, double *y) |
49 | { |
50 | union {double f; uint64_t i;} u = {x}; |
51 | double_t z,w,t,r,fn; |
52 | double tx[3],ty[2]; |
53 | uint32_t ix; |
54 | int sign, n, ex, ey, i; |
55 | |
56 | sign = u.i>>63; |
57 | ix = u.i>>32 & 0x7fffffff; |
58 | if (ix <= 0x400f6a7a) { |
| 1 | Assuming 'ix' is > 1074752122 | |
|
| |
59 | if ((ix & 0xfffff) == 0x921fb) |
60 | goto medium; |
61 | if (ix <= 0x4002d97c) { |
62 | if (!sign) { |
63 | z = x - pio2_1; |
64 | y[0] = z - pio2_1t; |
65 | y[1] = (z-y[0]) - pio2_1t; |
66 | return 1; |
67 | } else { |
68 | z = x + pio2_1; |
69 | y[0] = z + pio2_1t; |
70 | y[1] = (z-y[0]) + pio2_1t; |
71 | return -1; |
72 | } |
73 | } else { |
74 | if (!sign) { |
75 | z = x - 2*pio2_1; |
76 | y[0] = z - 2*pio2_1t; |
77 | y[1] = (z-y[0]) - 2*pio2_1t; |
78 | return 2; |
79 | } else { |
80 | z = x + 2*pio2_1; |
81 | y[0] = z + 2*pio2_1t; |
82 | y[1] = (z-y[0]) + 2*pio2_1t; |
83 | return -2; |
84 | } |
85 | } |
86 | } |
87 | if (ix <= 0x401c463b) { |
| 3 | | Assuming 'ix' is > 1075594811 | |
|
| |
88 | if (ix <= 0x4015fdbc) { |
89 | if (ix == 0x4012d97c) |
90 | goto medium; |
91 | if (!sign) { |
92 | z = x - 3*pio2_1; |
93 | y[0] = z - 3*pio2_1t; |
94 | y[1] = (z-y[0]) - 3*pio2_1t; |
95 | return 3; |
96 | } else { |
97 | z = x + 3*pio2_1; |
98 | y[0] = z + 3*pio2_1t; |
99 | y[1] = (z-y[0]) + 3*pio2_1t; |
100 | return -3; |
101 | } |
102 | } else { |
103 | if (ix == 0x401921fb) |
104 | goto medium; |
105 | if (!sign) { |
106 | z = x - 4*pio2_1; |
107 | y[0] = z - 4*pio2_1t; |
108 | y[1] = (z-y[0]) - 4*pio2_1t; |
109 | return 4; |
110 | } else { |
111 | z = x + 4*pio2_1; |
112 | y[0] = z + 4*pio2_1t; |
113 | y[1] = (z-y[0]) + 4*pio2_1t; |
114 | return -4; |
115 | } |
116 | } |
117 | } |
118 | if (ix < 0x413921fb) { |
| 5 | | Assuming 'ix' is >= 1094263291 | |
|
| |
119 | medium: |
120 | |
121 | fn = (double_t)x*invpio2 + toint - toint; |
122 | n = (int32_t)fn; |
123 | r = x - fn*pio2_1; |
124 | w = fn*pio2_1t; |
125 | y[0] = r - w; |
126 | u.f = y[0]; |
127 | ey = u.i>>52 & 0x7ff; |
128 | ex = ix>>20; |
129 | if (ex - ey > 16) { |
130 | t = r; |
131 | w = fn*pio2_2; |
132 | r = t - w; |
133 | w = fn*pio2_2t - ((t-r)-w); |
134 | y[0] = r - w; |
135 | u.f = y[0]; |
136 | ey = u.i>>52 & 0x7ff; |
137 | if (ex - ey > 49) { |
138 | t = r; |
139 | w = fn*pio2_3; |
140 | r = t - w; |
141 | w = fn*pio2_3t - ((t-r)-w); |
142 | y[0] = r - w; |
143 | } |
144 | } |
145 | y[1] = (r - y[0]) - w; |
146 | return n; |
147 | } |
148 | |
149 | |
150 | |
151 | if (ix >= 0x7ff00000) { |
| 7 | | Assuming 'ix' is < 2146435072 | |
|
| |
152 | y[0] = y[1] = x - x; |
153 | return 0; |
154 | } |
155 | |
156 | u.f = x; |
157 | u.i &= (uint64_t)-1>>12; |
158 | u.i |= (uint64_t)(0x3ff + 23)<<52; |
159 | z = u.f; |
160 | for (i=0; i < 2; i++) { |
| 9 | | Loop condition is true. Entering loop body | |
|
| 10 | | Loop condition is true. Entering loop body | |
|
| 11 | | Loop condition is false. Execution continues on line 164 | |
|
161 | tx[i] = (double)(int32_t)z; |
162 | z = (z-tx[i])*0x1p24; |
163 | } |
164 | tx[i] = z; |
165 | |
166 | while (tx[i] == 0.0) |
| 12 | | Loop condition is true. Entering loop body | |
|
| 13 | | Loop condition is true. Entering loop body | |
|
| 14 | | Loop condition is true. Entering loop body | |
|
| 15 | | The left operand of '==' is a garbage value |
|
167 | i--; |
168 | n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1); |
169 | if (sign) { |
170 | y[0] = -ty[0]; |
171 | y[1] = -ty[1]; |
172 | return -n; |
173 | } |
174 | y[0] = ty[0]; |
175 | y[1] = ty[1]; |
176 | return n; |
177 | } |