/src/moddable/xs/tools/fdlibm/e_rem_pio2.c
Line | Count | Source |
1 | | |
2 | | /* |
3 | | * ==================================================== |
4 | | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 | | * |
6 | | * Developed at SunSoft, a Sun Microsystems, Inc. business. |
7 | | * Permission to use, copy, modify, and distribute this |
8 | | * software is freely granted, provided that this notice |
9 | | * is preserved. |
10 | | * ==================================================== |
11 | | * |
12 | | * Optimized by Bruce D. Evans. |
13 | | */ |
14 | | |
15 | | /* __ieee754_rem_pio2(x,y) |
16 | | * |
17 | | * return the remainder of x rem pi/2 in y[0]+y[1] |
18 | | * use __kernel_rem_pio2() |
19 | | */ |
20 | | |
21 | | #include "math_private.h" |
22 | | |
23 | | /* |
24 | | * invpio2: 53 bits of 2/pi |
25 | | * pio2_1: first 33 bit of pi/2 |
26 | | * pio2_1t: pi/2 - pio2_1 |
27 | | * pio2_2: second 33 bit of pi/2 |
28 | | * pio2_2t: pi/2 - (pio2_1+pio2_2) |
29 | | * pio2_3: third 33 bit of pi/2 |
30 | | * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) |
31 | | */ |
32 | | |
33 | | static const double |
34 | | zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ |
35 | | two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ |
36 | | invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ |
37 | | pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ |
38 | | pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ |
39 | | pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ |
40 | | pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ |
41 | | pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ |
42 | | pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ |
43 | | |
44 | | #ifdef INLINE_REM_PIO2 |
45 | | static inline |
46 | | #endif |
47 | | int |
48 | | __ieee754_rem_pio2(double x, double *y) |
49 | 108k | { |
50 | 108k | double z,w,t,r,fn; |
51 | 108k | double tx[3],ty[2]; |
52 | 108k | int32_t e0,i,j,nx,n,ix,hx; |
53 | 108k | u_int32_t low; |
54 | | |
55 | 108k | GET_HIGH_WORD(hx,x); /* high word of x */ |
56 | 108k | ix = hx&0x7fffffff; |
57 | | #if 0 /* Must be handled in caller. */ |
58 | | if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ |
59 | | {y[0] = x; y[1] = 0; return 0;} |
60 | | #endif |
61 | 108k | if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */ |
62 | 20.3k | if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */ |
63 | 0 | goto medium; /* cancellation -- use medium case */ |
64 | 20.3k | if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */ |
65 | 12.6k | if (hx > 0) { |
66 | 10.7k | z = x - pio2_1; /* one round good to 85 bits */ |
67 | 10.7k | y[0] = z - pio2_1t; |
68 | 10.7k | y[1] = (z-y[0])-pio2_1t; |
69 | 10.7k | return 1; |
70 | 10.7k | } else { |
71 | 1.86k | z = x + pio2_1; |
72 | 1.86k | y[0] = z + pio2_1t; |
73 | 1.86k | y[1] = (z-y[0])+pio2_1t; |
74 | 1.86k | return -1; |
75 | 1.86k | } |
76 | 12.6k | } else { |
77 | 7.79k | if (hx > 0) { |
78 | 5.22k | z = x - 2*pio2_1; |
79 | 5.22k | y[0] = z - 2*pio2_1t; |
80 | 5.22k | y[1] = (z-y[0])-2*pio2_1t; |
81 | 5.22k | return 2; |
82 | 5.22k | } else { |
83 | 2.56k | z = x + 2*pio2_1; |
84 | 2.56k | y[0] = z + 2*pio2_1t; |
85 | 2.56k | y[1] = (z-y[0])+2*pio2_1t; |
86 | 2.56k | return -2; |
87 | 2.56k | } |
88 | 7.79k | } |
89 | 20.3k | } |
90 | 87.8k | if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */ |
91 | 8.87k | if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */ |
92 | 5.92k | if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */ |
93 | 0 | goto medium; |
94 | 5.92k | if (hx > 0) { |
95 | 3.97k | z = x - 3*pio2_1; |
96 | 3.97k | y[0] = z - 3*pio2_1t; |
97 | 3.97k | y[1] = (z-y[0])-3*pio2_1t; |
98 | 3.97k | return 3; |
99 | 3.97k | } else { |
100 | 1.94k | z = x + 3*pio2_1; |
101 | 1.94k | y[0] = z + 3*pio2_1t; |
102 | 1.94k | y[1] = (z-y[0])+3*pio2_1t; |
103 | 1.94k | return -3; |
104 | 1.94k | } |
105 | 5.92k | } else { |
106 | 2.94k | if (ix == 0x401921fb) /* |x| ~= 4pi/2 */ |
107 | 0 | goto medium; |
108 | 2.94k | if (hx > 0) { |
109 | 1.70k | z = x - 4*pio2_1; |
110 | 1.70k | y[0] = z - 4*pio2_1t; |
111 | 1.70k | y[1] = (z-y[0])-4*pio2_1t; |
112 | 1.70k | return 4; |
113 | 1.70k | } else { |
114 | 1.24k | z = x + 4*pio2_1; |
115 | 1.24k | y[0] = z + 4*pio2_1t; |
116 | 1.24k | y[1] = (z-y[0])+4*pio2_1t; |
117 | 1.24k | return -4; |
118 | 1.24k | } |
119 | 2.94k | } |
120 | 8.87k | } |
121 | 78.9k | if(ix<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ |
122 | 71.2k | medium: |
123 | 71.2k | fn = rnint((double_t)x*invpio2); |
124 | 71.2k | n = irint(fn); |
125 | 71.2k | r = x-fn*pio2_1; |
126 | 71.2k | w = fn*pio2_1t; /* 1st round good to 85 bit */ |
127 | 71.2k | { |
128 | 71.2k | u_int32_t high; |
129 | 71.2k | j = ix>>20; |
130 | 71.2k | y[0] = r-w; |
131 | 71.2k | GET_HIGH_WORD(high,y[0]); |
132 | 71.2k | i = j-((high>>20)&0x7ff); |
133 | 71.2k | if(i>16) { /* 2nd iteration needed, good to 118 */ |
134 | 4.80k | t = r; |
135 | 4.80k | w = fn*pio2_2; |
136 | 4.80k | r = t-w; |
137 | 4.80k | w = fn*pio2_2t-((t-r)-w); |
138 | 4.80k | y[0] = r-w; |
139 | 4.80k | GET_HIGH_WORD(high,y[0]); |
140 | 4.80k | i = j-((high>>20)&0x7ff); |
141 | 4.80k | if(i>49) { /* 3rd iteration need, 151 bits acc */ |
142 | 0 | t = r; /* will cover all possible cases */ |
143 | 0 | w = fn*pio2_3; |
144 | 0 | r = t-w; |
145 | 0 | w = fn*pio2_3t-((t-r)-w); |
146 | 0 | y[0] = r-w; |
147 | 0 | } |
148 | 4.80k | } |
149 | 71.2k | } |
150 | 71.2k | y[1] = (r-y[0])-w; |
151 | 71.2k | return n; |
152 | 71.2k | } |
153 | | /* |
154 | | * all other (large) arguments |
155 | | */ |
156 | 7.77k | if(ix>=0x7ff00000) { /* x is inf or NaN */ |
157 | 0 | y[0]=y[1]=x-x; return 0; |
158 | 0 | } |
159 | | /* set z = scalbn(|x|,ilogb(x)-23) */ |
160 | 7.77k | GET_LOW_WORD(low,x); |
161 | 7.77k | e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ |
162 | 7.77k | INSERT_WORDS(z, ix - ((int32_t)((u_int32_t)e0<<20)), low); |
163 | 23.3k | for(i=0;i<2;i++) { |
164 | 15.5k | tx[i] = (double)((int32_t)(z)); |
165 | 15.5k | z = (z-tx[i])*two24; |
166 | 15.5k | } |
167 | 7.77k | tx[2] = z; |
168 | 7.77k | nx = 3; |
169 | 11.4k | while(tx[nx-1]==zero) nx--; /* skip zero term */ |
170 | 7.77k | n = __kernel_rem_pio2(tx,ty,e0,nx,1); |
171 | 7.77k | if(hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;} |
172 | 6.17k | y[0] = ty[0]; y[1] = ty[1]; return n; |
173 | 7.77k | } Unexecuted instantiation: __ieee754_rem_pio2 s_cos.c:__ieee754_rem_pio2 Line | Count | Source | 49 | 24.1k | { | 50 | 24.1k | double z,w,t,r,fn; | 51 | 24.1k | double tx[3],ty[2]; | 52 | 24.1k | int32_t e0,i,j,nx,n,ix,hx; | 53 | 24.1k | u_int32_t low; | 54 | | | 55 | 24.1k | GET_HIGH_WORD(hx,x); /* high word of x */ | 56 | 24.1k | ix = hx&0x7fffffff; | 57 | | #if 0 /* Must be handled in caller. */ | 58 | | if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ | 59 | | {y[0] = x; y[1] = 0; return 0;} | 60 | | #endif | 61 | 24.1k | if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */ | 62 | 11.2k | if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */ | 63 | 0 | goto medium; /* cancellation -- use medium case */ | 64 | 11.2k | if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */ | 65 | 9.32k | if (hx > 0) { | 66 | 8.08k | z = x - pio2_1; /* one round good to 85 bits */ | 67 | 8.08k | y[0] = z - pio2_1t; | 68 | 8.08k | y[1] = (z-y[0])-pio2_1t; | 69 | 8.08k | return 1; | 70 | 8.08k | } else { | 71 | 1.23k | z = x + pio2_1; | 72 | 1.23k | y[0] = z + pio2_1t; | 73 | 1.23k | y[1] = (z-y[0])+pio2_1t; | 74 | 1.23k | return -1; | 75 | 1.23k | } | 76 | 9.32k | } else { | 77 | 1.92k | if (hx > 0) { | 78 | 1.06k | z = x - 2*pio2_1; | 79 | 1.06k | y[0] = z - 2*pio2_1t; | 80 | 1.06k | y[1] = (z-y[0])-2*pio2_1t; | 81 | 1.06k | return 2; | 82 | 1.06k | } else { | 83 | 859 | z = x + 2*pio2_1; | 84 | 859 | y[0] = z + 2*pio2_1t; | 85 | 859 | y[1] = (z-y[0])+2*pio2_1t; | 86 | 859 | return -2; | 87 | 859 | } | 88 | 1.92k | } | 89 | 11.2k | } | 90 | 12.8k | if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */ | 91 | 3.13k | if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */ | 92 | 1.80k | if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */ | 93 | 0 | goto medium; | 94 | 1.80k | if (hx > 0) { | 95 | 1.03k | z = x - 3*pio2_1; | 96 | 1.03k | y[0] = z - 3*pio2_1t; | 97 | 1.03k | y[1] = (z-y[0])-3*pio2_1t; | 98 | 1.03k | return 3; | 99 | 1.03k | } else { | 100 | 776 | z = x + 3*pio2_1; | 101 | 776 | y[0] = z + 3*pio2_1t; | 102 | 776 | y[1] = (z-y[0])+3*pio2_1t; | 103 | 776 | return -3; | 104 | 776 | } | 105 | 1.80k | } else { | 106 | 1.32k | if (ix == 0x401921fb) /* |x| ~= 4pi/2 */ | 107 | 0 | goto medium; | 108 | 1.32k | if (hx > 0) { | 109 | 707 | z = x - 4*pio2_1; | 110 | 707 | y[0] = z - 4*pio2_1t; | 111 | 707 | y[1] = (z-y[0])-4*pio2_1t; | 112 | 707 | return 4; | 113 | 707 | } else { | 114 | 618 | z = x + 4*pio2_1; | 115 | 618 | y[0] = z + 4*pio2_1t; | 116 | 618 | y[1] = (z-y[0])+4*pio2_1t; | 117 | 618 | return -4; | 118 | 618 | } | 119 | 1.32k | } | 120 | 3.13k | } | 121 | 9.73k | if(ix<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ | 122 | 4.66k | medium: | 123 | 4.66k | fn = rnint((double_t)x*invpio2); | 124 | 4.66k | n = irint(fn); | 125 | 4.66k | r = x-fn*pio2_1; | 126 | 4.66k | w = fn*pio2_1t; /* 1st round good to 85 bit */ | 127 | 4.66k | { | 128 | 4.66k | u_int32_t high; | 129 | 4.66k | j = ix>>20; | 130 | 4.66k | y[0] = r-w; | 131 | 4.66k | GET_HIGH_WORD(high,y[0]); | 132 | 4.66k | i = j-((high>>20)&0x7ff); | 133 | 4.66k | if(i>16) { /* 2nd iteration needed, good to 118 */ | 134 | 561 | t = r; | 135 | 561 | w = fn*pio2_2; | 136 | 561 | r = t-w; | 137 | 561 | w = fn*pio2_2t-((t-r)-w); | 138 | 561 | y[0] = r-w; | 139 | 561 | GET_HIGH_WORD(high,y[0]); | 140 | 561 | i = j-((high>>20)&0x7ff); | 141 | 561 | if(i>49) { /* 3rd iteration need, 151 bits acc */ | 142 | 0 | t = r; /* will cover all possible cases */ | 143 | 0 | w = fn*pio2_3; | 144 | 0 | r = t-w; | 145 | 0 | w = fn*pio2_3t-((t-r)-w); | 146 | 0 | y[0] = r-w; | 147 | 0 | } | 148 | 561 | } | 149 | 4.66k | } | 150 | 4.66k | y[1] = (r-y[0])-w; | 151 | 4.66k | return n; | 152 | 4.66k | } | 153 | | /* | 154 | | * all other (large) arguments | 155 | | */ | 156 | 5.06k | if(ix>=0x7ff00000) { /* x is inf or NaN */ | 157 | 0 | y[0]=y[1]=x-x; return 0; | 158 | 0 | } | 159 | | /* set z = scalbn(|x|,ilogb(x)-23) */ | 160 | 5.06k | GET_LOW_WORD(low,x); | 161 | 5.06k | e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ | 162 | 5.06k | INSERT_WORDS(z, ix - ((int32_t)((u_int32_t)e0<<20)), low); | 163 | 15.2k | for(i=0;i<2;i++) { | 164 | 10.1k | tx[i] = (double)((int32_t)(z)); | 165 | 10.1k | z = (z-tx[i])*two24; | 166 | 10.1k | } | 167 | 5.06k | tx[2] = z; | 168 | 5.06k | nx = 3; | 169 | 6.03k | while(tx[nx-1]==zero) nx--; /* skip zero term */ | 170 | 5.06k | n = __kernel_rem_pio2(tx,ty,e0,nx,1); | 171 | 5.06k | if(hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;} | 172 | 4.26k | y[0] = ty[0]; y[1] = ty[1]; return n; | 173 | 5.06k | } |
s_sin.c:__ieee754_rem_pio2 Line | Count | Source | 49 | 8.74k | { | 50 | 8.74k | double z,w,t,r,fn; | 51 | 8.74k | double tx[3],ty[2]; | 52 | 8.74k | int32_t e0,i,j,nx,n,ix,hx; | 53 | 8.74k | u_int32_t low; | 54 | | | 55 | 8.74k | GET_HIGH_WORD(hx,x); /* high word of x */ | 56 | 8.74k | ix = hx&0x7fffffff; | 57 | | #if 0 /* Must be handled in caller. */ | 58 | | if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ | 59 | | {y[0] = x; y[1] = 0; return 0;} | 60 | | #endif | 61 | 8.74k | if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */ | 62 | 3.67k | if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */ | 63 | 0 | goto medium; /* cancellation -- use medium case */ | 64 | 3.67k | if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */ | 65 | 1.40k | if (hx > 0) { | 66 | 985 | z = x - pio2_1; /* one round good to 85 bits */ | 67 | 985 | y[0] = z - pio2_1t; | 68 | 985 | y[1] = (z-y[0])-pio2_1t; | 69 | 985 | return 1; | 70 | 985 | } else { | 71 | 423 | z = x + pio2_1; | 72 | 423 | y[0] = z + pio2_1t; | 73 | 423 | y[1] = (z-y[0])+pio2_1t; | 74 | 423 | return -1; | 75 | 423 | } | 76 | 2.26k | } else { | 77 | 2.26k | if (hx > 0) { | 78 | 783 | z = x - 2*pio2_1; | 79 | 783 | y[0] = z - 2*pio2_1t; | 80 | 783 | y[1] = (z-y[0])-2*pio2_1t; | 81 | 783 | return 2; | 82 | 1.48k | } else { | 83 | 1.48k | z = x + 2*pio2_1; | 84 | 1.48k | y[0] = z + 2*pio2_1t; | 85 | 1.48k | y[1] = (z-y[0])+2*pio2_1t; | 86 | 1.48k | return -2; | 87 | 1.48k | } | 88 | 2.26k | } | 89 | 3.67k | } | 90 | 5.06k | if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */ | 91 | 2.63k | if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */ | 92 | 1.88k | if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */ | 93 | 0 | goto medium; | 94 | 1.88k | if (hx > 0) { | 95 | 921 | z = x - 3*pio2_1; | 96 | 921 | y[0] = z - 3*pio2_1t; | 97 | 921 | y[1] = (z-y[0])-3*pio2_1t; | 98 | 921 | return 3; | 99 | 959 | } else { | 100 | 959 | z = x + 3*pio2_1; | 101 | 959 | y[0] = z + 3*pio2_1t; | 102 | 959 | y[1] = (z-y[0])+3*pio2_1t; | 103 | 959 | return -3; | 104 | 959 | } | 105 | 1.88k | } else { | 106 | 759 | if (ix == 0x401921fb) /* |x| ~= 4pi/2 */ | 107 | 0 | goto medium; | 108 | 759 | if (hx > 0) { | 109 | 336 | z = x - 4*pio2_1; | 110 | 336 | y[0] = z - 4*pio2_1t; | 111 | 336 | y[1] = (z-y[0])-4*pio2_1t; | 112 | 336 | return 4; | 113 | 423 | } else { | 114 | 423 | z = x + 4*pio2_1; | 115 | 423 | y[0] = z + 4*pio2_1t; | 116 | 423 | y[1] = (z-y[0])+4*pio2_1t; | 117 | 423 | return -4; | 118 | 423 | } | 119 | 759 | } | 120 | 2.63k | } | 121 | 2.43k | if(ix<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ | 122 | 979 | medium: | 123 | 979 | fn = rnint((double_t)x*invpio2); | 124 | 979 | n = irint(fn); | 125 | 979 | r = x-fn*pio2_1; | 126 | 979 | w = fn*pio2_1t; /* 1st round good to 85 bit */ | 127 | 979 | { | 128 | 979 | u_int32_t high; | 129 | 979 | j = ix>>20; | 130 | 979 | y[0] = r-w; | 131 | 979 | GET_HIGH_WORD(high,y[0]); | 132 | 979 | i = j-((high>>20)&0x7ff); | 133 | 979 | if(i>16) { /* 2nd iteration needed, good to 118 */ | 134 | 239 | t = r; | 135 | 239 | w = fn*pio2_2; | 136 | 239 | r = t-w; | 137 | 239 | w = fn*pio2_2t-((t-r)-w); | 138 | 239 | y[0] = r-w; | 139 | 239 | GET_HIGH_WORD(high,y[0]); | 140 | 239 | i = j-((high>>20)&0x7ff); | 141 | 239 | if(i>49) { /* 3rd iteration need, 151 bits acc */ | 142 | 0 | t = r; /* will cover all possible cases */ | 143 | 0 | w = fn*pio2_3; | 144 | 0 | r = t-w; | 145 | 0 | w = fn*pio2_3t-((t-r)-w); | 146 | 0 | y[0] = r-w; | 147 | 0 | } | 148 | 239 | } | 149 | 979 | } | 150 | 979 | y[1] = (r-y[0])-w; | 151 | 979 | return n; | 152 | 979 | } | 153 | | /* | 154 | | * all other (large) arguments | 155 | | */ | 156 | 1.45k | if(ix>=0x7ff00000) { /* x is inf or NaN */ | 157 | 0 | y[0]=y[1]=x-x; return 0; | 158 | 0 | } | 159 | | /* set z = scalbn(|x|,ilogb(x)-23) */ | 160 | 1.45k | GET_LOW_WORD(low,x); | 161 | 1.45k | e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ | 162 | 1.45k | INSERT_WORDS(z, ix - ((int32_t)((u_int32_t)e0<<20)), low); | 163 | 4.35k | for(i=0;i<2;i++) { | 164 | 2.90k | tx[i] = (double)((int32_t)(z)); | 165 | 2.90k | z = (z-tx[i])*two24; | 166 | 2.90k | } | 167 | 1.45k | tx[2] = z; | 168 | 1.45k | nx = 3; | 169 | 3.50k | while(tx[nx-1]==zero) nx--; /* skip zero term */ | 170 | 1.45k | n = __kernel_rem_pio2(tx,ty,e0,nx,1); | 171 | 1.45k | if(hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;} | 172 | 1.18k | y[0] = ty[0]; y[1] = ty[1]; return n; | 173 | 1.45k | } |
s_tan.c:__ieee754_rem_pio2 Line | Count | Source | 49 | 75.4k | { | 50 | 75.4k | double z,w,t,r,fn; | 51 | 75.4k | double tx[3],ty[2]; | 52 | 75.4k | int32_t e0,i,j,nx,n,ix,hx; | 53 | 75.4k | u_int32_t low; | 54 | | | 55 | 75.4k | GET_HIGH_WORD(hx,x); /* high word of x */ | 56 | 75.4k | ix = hx&0x7fffffff; | 57 | | #if 0 /* Must be handled in caller. */ | 58 | | if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ | 59 | | {y[0] = x; y[1] = 0; return 0;} | 60 | | #endif | 61 | 75.4k | if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */ | 62 | 5.48k | if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */ | 63 | 0 | goto medium; /* cancellation -- use medium case */ | 64 | 5.48k | if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */ | 65 | 1.87k | if (hx > 0) { | 66 | 1.66k | z = x - pio2_1; /* one round good to 85 bits */ | 67 | 1.66k | y[0] = z - pio2_1t; | 68 | 1.66k | y[1] = (z-y[0])-pio2_1t; | 69 | 1.66k | return 1; | 70 | 1.66k | } else { | 71 | 206 | z = x + pio2_1; | 72 | 206 | y[0] = z + pio2_1t; | 73 | 206 | y[1] = (z-y[0])+pio2_1t; | 74 | 206 | return -1; | 75 | 206 | } | 76 | 3.60k | } else { | 77 | 3.60k | if (hx > 0) { | 78 | 3.38k | z = x - 2*pio2_1; | 79 | 3.38k | y[0] = z - 2*pio2_1t; | 80 | 3.38k | y[1] = (z-y[0])-2*pio2_1t; | 81 | 3.38k | return 2; | 82 | 3.38k | } else { | 83 | 225 | z = x + 2*pio2_1; | 84 | 225 | y[0] = z + 2*pio2_1t; | 85 | 225 | y[1] = (z-y[0])+2*pio2_1t; | 86 | 225 | return -2; | 87 | 225 | } | 88 | 3.60k | } | 89 | 5.48k | } | 90 | 69.9k | if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */ | 91 | 3.10k | if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */ | 92 | 2.23k | if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */ | 93 | 0 | goto medium; | 94 | 2.23k | if (hx > 0) { | 95 | 2.02k | z = x - 3*pio2_1; | 96 | 2.02k | y[0] = z - 3*pio2_1t; | 97 | 2.02k | y[1] = (z-y[0])-3*pio2_1t; | 98 | 2.02k | return 3; | 99 | 2.02k | } else { | 100 | 210 | z = x + 3*pio2_1; | 101 | 210 | y[0] = z + 3*pio2_1t; | 102 | 210 | y[1] = (z-y[0])+3*pio2_1t; | 103 | 210 | return -3; | 104 | 210 | } | 105 | 2.23k | } else { | 106 | 862 | if (ix == 0x401921fb) /* |x| ~= 4pi/2 */ | 107 | 0 | goto medium; | 108 | 862 | if (hx > 0) { | 109 | 657 | z = x - 4*pio2_1; | 110 | 657 | y[0] = z - 4*pio2_1t; | 111 | 657 | y[1] = (z-y[0])-4*pio2_1t; | 112 | 657 | return 4; | 113 | 657 | } else { | 114 | 205 | z = x + 4*pio2_1; | 115 | 205 | y[0] = z + 4*pio2_1t; | 116 | 205 | y[1] = (z-y[0])+4*pio2_1t; | 117 | 205 | return -4; | 118 | 205 | } | 119 | 862 | } | 120 | 3.10k | } | 121 | 66.8k | if(ix<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ | 122 | 65.5k | medium: | 123 | 65.5k | fn = rnint((double_t)x*invpio2); | 124 | 65.5k | n = irint(fn); | 125 | 65.5k | r = x-fn*pio2_1; | 126 | 65.5k | w = fn*pio2_1t; /* 1st round good to 85 bit */ | 127 | 65.5k | { | 128 | 65.5k | u_int32_t high; | 129 | 65.5k | j = ix>>20; | 130 | 65.5k | y[0] = r-w; | 131 | 65.5k | GET_HIGH_WORD(high,y[0]); | 132 | 65.5k | i = j-((high>>20)&0x7ff); | 133 | 65.5k | if(i>16) { /* 2nd iteration needed, good to 118 */ | 134 | 4.00k | t = r; | 135 | 4.00k | w = fn*pio2_2; | 136 | 4.00k | r = t-w; | 137 | 4.00k | w = fn*pio2_2t-((t-r)-w); | 138 | 4.00k | y[0] = r-w; | 139 | 4.00k | GET_HIGH_WORD(high,y[0]); | 140 | 4.00k | i = j-((high>>20)&0x7ff); | 141 | 4.00k | if(i>49) { /* 3rd iteration need, 151 bits acc */ | 142 | 0 | t = r; /* will cover all possible cases */ | 143 | 0 | w = fn*pio2_3; | 144 | 0 | r = t-w; | 145 | 0 | w = fn*pio2_3t-((t-r)-w); | 146 | 0 | y[0] = r-w; | 147 | 0 | } | 148 | 4.00k | } | 149 | 65.5k | } | 150 | 65.5k | y[1] = (r-y[0])-w; | 151 | 65.5k | return n; | 152 | 65.5k | } | 153 | | /* | 154 | | * all other (large) arguments | 155 | | */ | 156 | 1.25k | if(ix>=0x7ff00000) { /* x is inf or NaN */ | 157 | 0 | y[0]=y[1]=x-x; return 0; | 158 | 0 | } | 159 | | /* set z = scalbn(|x|,ilogb(x)-23) */ | 160 | 1.25k | GET_LOW_WORD(low,x); | 161 | 1.25k | e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ | 162 | 1.25k | INSERT_WORDS(z, ix - ((int32_t)((u_int32_t)e0<<20)), low); | 163 | 3.77k | for(i=0;i<2;i++) { | 164 | 2.51k | tx[i] = (double)((int32_t)(z)); | 165 | 2.51k | z = (z-tx[i])*two24; | 166 | 2.51k | } | 167 | 1.25k | tx[2] = z; | 168 | 1.25k | nx = 3; | 169 | 1.87k | while(tx[nx-1]==zero) nx--; /* skip zero term */ | 170 | 1.25k | n = __kernel_rem_pio2(tx,ty,e0,nx,1); | 171 | 1.25k | if(hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;} | 172 | 731 | y[0] = ty[0]; y[1] = ty[1]; return n; | 173 | 1.25k | } |
|