Mercurial > sdl-ios-xcode

comparison src/video/e_pow.h @ 1895:c121d94672cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
SDL 1.2 is moving to a branch, and SDL 1.3 is becoming the head.
author Sam Lantinga <slouken@libsdl.org>
date Mon, 10 Jul 2006 21:04:37 +0000
parents 7a610f25c12f
children edd2839b36f7
comparison
equal deleted inserted replaced
1894:c69cee13dd76 1895:c121d94672cb
65 #ifdef __STDC__ 65 #ifdef __STDC__
66 static const double 66 static const double
67 #else 67 #else
68 static double 68 static double
69 #endif 69 #endif
70 bp[] = {1.0, 1.5,}, 70 bp[] = { 1.0, 1.5, }, dp_h[] = {
71 dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ 71 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
72 dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ 72
73 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ 73 dp_l[] = {
74 L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ 74 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
75 L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ 75
76 L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ 76 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
77 L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ 77 L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
78 L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ 78 L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
79 L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ 79 L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
80 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ 80 L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
81 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ 81 L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
82 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ 82 L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
83 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ 83 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
84 P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */ 84 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
85 lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ 85 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
86 lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ 86 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
87 lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ 87 P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
88 ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ 88 lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
89 cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ 89 lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
90 cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ 90 lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
91 cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ 91 ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
92 ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ 92 cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
93 ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ 93 cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
94 ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ 94 cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h */
95 ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
96 ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2 */
97 ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail */
95 98
96 #ifdef __STDC__ 99 #ifdef __STDC__
97 double __ieee754_pow(double x, double y) 100 double
101 __ieee754_pow(double x, double y)
98 #else 102 #else
99 double __ieee754_pow(x,y) 103 double
100 double x, y; 104 __ieee754_pow(x, y)
105 double x, y;
101 #endif 106 #endif
102 { 107 {
103 double z,ax,z_h,z_l,p_h,p_l; 108 double z, ax, z_h, z_l, p_h, p_l;
104 double y1,t1,t2,r,s,t,u,v,w; 109 double y1, t1, t2, r, s, t, u, v, w;
105 int32_t i,j,k,yisint,n; 110 int32_t i, j, k, yisint, n;
106 int32_t hx,hy,ix,iy; 111 int32_t hx, hy, ix, iy;
107 u_int32_t lx,ly; 112 u_int32_t lx, ly;
108 113
109 EXTRACT_WORDS(hx,lx,x); 114 EXTRACT_WORDS(hx, lx, x);
110 EXTRACT_WORDS(hy,ly,y); 115 EXTRACT_WORDS(hy, ly, y);
111 ix = hx&0x7fffffff; iy = hy&0x7fffffff; 116 ix = hx & 0x7fffffff;
117 iy = hy & 0x7fffffff;
112 118
113 /* y==zero: x**0 = 1 */ 119 /* y==zero: x**0 = 1 */
114 if((iy|ly)==0) return one; 120 if ((iy | ly) == 0)
121 return one;
115 122
116 /* +-NaN return x+y */ 123 /* +-NaN return x+y */
117 if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || 124 if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) ||
118 iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) 125 iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0)))
119 return x+y; 126 return x + y;
120 127
121 /* determine if y is an odd int when x < 0 128 /* determine if y is an odd int when x < 0
122 * yisint = 0 ... y is not an integer 129 * yisint = 0 ... y is not an integer
123 * yisint = 1 ... y is an odd int 130 * yisint = 1 ... y is an odd int
124 * yisint = 2 ... y is an even int 131 * yisint = 2 ... y is an even int
125 */ 132 */
126 yisint = 0; 133 yisint = 0;
127 if(hx<0) { 134 if (hx < 0) {
128 if(iy>=0x43400000) yisint = 2; /* even integer y */ 135 if (iy >= 0x43400000)
129 else if(iy>=0x3ff00000) { 136 yisint = 2; /* even integer y */
130 k = (iy>>20)-0x3ff; /* exponent */ 137 else if (iy >= 0x3ff00000) {
131 if(k>20) { 138 k = (iy >> 20) - 0x3ff; /* exponent */
132 j = ly>>(52-k); 139 if (k > 20) {
133 if((u_int32_t)(j<<(52-k))==ly) yisint = 2-(j&1); 140 j = ly >> (52 - k);
134 } else if(ly==0) { 141 if ((u_int32_t) (j << (52 - k)) == ly)
135 j = iy>>(20-k); 142 yisint = 2 - (j & 1);
136 if((j<<(20-k))==iy) yisint = 2-(j&1); 143 } else if (ly == 0) {
137 } 144 j = iy >> (20 - k);
138 } 145 if ((j << (20 - k)) == iy)
139 } 146 yisint = 2 - (j & 1);
147 }
148 }
149 }
140 150
141 /* special value of y */ 151 /* special value of y */
142 if(ly==0) { 152 if (ly == 0) {
143 if (iy==0x7ff00000) { /* y is +-inf */ 153 if (iy == 0x7ff00000) { /* y is +-inf */
144 if(((ix-0x3ff00000)|lx)==0) 154 if (((ix - 0x3ff00000) | lx) == 0)
145 return y - y; /* inf**+-1 is NaN */ 155 return y - y; /* inf**+-1 is NaN */
146 else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */ 156 else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
147 return (hy>=0)? y: zero; 157 return (hy >= 0) ? y : zero;
148 else /* (|x|<1)**-,+inf = inf,0 */ 158 else /* (|x|<1)**-,+inf = inf,0 */
149 return (hy<0)?-y: zero; 159 return (hy < 0) ? -y : zero;
150 } 160 }
151 if(iy==0x3ff00000) { /* y is +-1 */ 161 if (iy == 0x3ff00000) { /* y is +-1 */
152 if(hy<0) return one/x; else return x; 162 if (hy < 0)
153 } 163 return one / x;
154 if(hy==0x40000000) return x*x; /* y is 2 */ 164 else
155 if(hy==0x3fe00000) { /* y is 0.5 */ 165 return x;
156 if(hx>=0) /* x >= +0 */ 166 }
157 return __ieee754_sqrt(x); 167 if (hy == 0x40000000)
158 } 168 return x * x; /* y is 2 */
159 } 169 if (hy == 0x3fe00000) { /* y is 0.5 */
160 170 if (hx >= 0) /* x >= +0 */
161 ax = x < 0 ? -x : x; /*fabs(x);*/ 171 return __ieee754_sqrt(x);
172 }
173 }
174
175 ax = x < 0 ? -x : x; /*fabs(x); */
162 /* special value of x */ 176 /* special value of x */
163 if(lx==0) { 177 if (lx == 0) {
164 if(ix==0x7ff00000||ix==0||ix==0x3ff00000){ 178 if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) {
165 z = ax; /*x is +-0,+-inf,+-1*/ 179 z = ax; /*x is +-0,+-inf,+-1 */
166 if(hy<0) z = one/z; /* z = (1/|x|) */ 180 if (hy < 0)
167 if(hx<0) { 181 z = one / z; /* z = (1/|x|) */
168 if(((ix-0x3ff00000)|yisint)==0) { 182 if (hx < 0) {
169 z = (z-z)/(z-z); /* (-1)**non-int is NaN */ 183 if (((ix - 0x3ff00000) | yisint) == 0) {
170 } else if(yisint==1) 184 z = (z - z) / (z - z); /* (-1)**non-int is NaN */
171 z = -z; /* (x<0)**odd = -(|x|**odd) */ 185 } else if (yisint == 1)
172 } 186 z = -z; /* (x<0)**odd = -(|x|**odd) */
173 return z; 187 }
174 } 188 return z;
175 } 189 }
190 }
176 191
177 /* (x<0)**(non-int) is NaN */ 192 /* (x<0)**(non-int) is NaN */
178 if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x); 193 if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
194 return (x - x) / (x - x);
179 195
180 /* |y| is huge */ 196 /* |y| is huge */
181 if(iy>0x41e00000) { /* if |y| > 2**31 */ 197 if (iy > 0x41e00000) { /* if |y| > 2**31 */
182 if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */ 198 if (iy > 0x43f00000) { /* if |y| > 2**64, must o/uflow */
183 if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny; 199 if (ix <= 0x3fefffff)
184 if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny; 200 return (hy < 0) ? huge * huge : tiny * tiny;
185 } 201 if (ix >= 0x3ff00000)
186 /* over/underflow if x is not close to one */ 202 return (hy > 0) ? huge * huge : tiny * tiny;
187 if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny; 203 }
188 if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny; 204 /* over/underflow if x is not close to one */
189 /* now |1-x| is tiny <= 2**-20, suffice to compute 205 if (ix < 0x3fefffff)
190 log(x) by x-x^2/2+x^3/3-x^4/4 */ 206 return (hy < 0) ? huge * huge : tiny * tiny;
191 t = x-1; /* t has 20 trailing zeros */ 207 if (ix > 0x3ff00000)
192 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25)); 208 return (hy > 0) ? huge * huge : tiny * tiny;
193 u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ 209 /* now |1-x| is tiny <= 2**-20, suffice to compute
194 v = t*ivln2_l-w*ivln2; 210 log(x) by x-x^2/2+x^3/3-x^4/4 */
195 t1 = u+v; 211 t = x - 1; /* t has 20 trailing zeros */
196 SET_LOW_WORD(t1,0); 212 w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25));
197 t2 = v-(t1-u); 213 u = ivln2_h * t; /* ivln2_h has 21 sig. bits */
198 } else { 214 v = t * ivln2_l - w * ivln2;
199 double s2,s_h,s_l,t_h,t_l; 215 t1 = u + v;
200 n = 0; 216 SET_LOW_WORD(t1, 0);
201 /* take care subnormal number */ 217 t2 = v - (t1 - u);
202 if(ix<0x00100000) 218 } else {
203 {ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); } 219 double s2, s_h, s_l, t_h, t_l;
204 n += ((ix)>>20)-0x3ff; 220 n = 0;
205 j = ix&0x000fffff; 221 /* take care subnormal number */
206 /* determine interval */ 222 if (ix < 0x00100000) {
207 ix = j|0x3ff00000; /* normalize ix */ 223 ax *= two53;
208 if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ 224 n -= 53;
209 else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ 225 GET_HIGH_WORD(ix, ax);
210 else {k=0;n+=1;ix -= 0x00100000;} 226 }
211 SET_HIGH_WORD(ax,ix); 227 n += ((ix) >> 20) - 0x3ff;
212 228 j = ix & 0x000fffff;
213 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ 229 /* determine interval */
214 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ 230 ix = j | 0x3ff00000; /* normalize ix */
215 v = one/(ax+bp[k]); 231 if (j <= 0x3988E)
216 s = u*v; 232 k = 0; /* |x|<sqrt(3/2) */
217 s_h = s; 233 else if (j < 0xBB67A)
218 SET_LOW_WORD(s_h,0); 234 k = 1; /* |x|<sqrt(3) */
219 /* t_h=ax+bp[k] High */ 235 else {
220 t_h = zero; 236 k = 0;
221 SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18)); 237 n += 1;
222 t_l = ax - (t_h-bp[k]); 238 ix -= 0x00100000;
223 s_l = v*((u-s_h*t_h)-s_h*t_l); 239 }
224 /* compute log(ax) */ 240 SET_HIGH_WORD(ax, ix);
225 s2 = s*s; 241
226 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); 242 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
227 r += s_l*(s_h+s); 243 u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
228 s2 = s_h*s_h; 244 v = one / (ax + bp[k]);
229 t_h = 3.0+s2+r; 245 s = u * v;
230 SET_LOW_WORD(t_h,0); 246 s_h = s;
231 t_l = r-((t_h-3.0)-s2); 247 SET_LOW_WORD(s_h, 0);
232 /* u+v = s*(1+...) */ 248 /* t_h=ax+bp[k] High */
233 u = s_h*t_h; 249 t_h = zero;
234 v = s_l*t_h+t_l*s; 250 SET_HIGH_WORD(t_h, ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18));
235 /* 2/(3log2)*(s+...) */ 251 t_l = ax - (t_h - bp[k]);
236 p_h = u+v; 252 s_l = v * ((u - s_h * t_h) - s_h * t_l);
237 SET_LOW_WORD(p_h,0); 253 /* compute log(ax) */
238 p_l = v-(p_h-u); 254 s2 = s * s;
239 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ 255 r = s2 * s2 * (L1 +
240 z_l = cp_l*p_h+p_l*cp+dp_l[k]; 256 s2 * (L2 +
241 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ 257 s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
242 t = (double)n; 258 r += s_l * (s_h + s);
243 t1 = (((z_h+z_l)+dp_h[k])+t); 259 s2 = s_h * s_h;
244 SET_LOW_WORD(t1,0); 260 t_h = 3.0 + s2 + r;
245 t2 = z_l-(((t1-t)-dp_h[k])-z_h); 261 SET_LOW_WORD(t_h, 0);
246 } 262 t_l = r - ((t_h - 3.0) - s2);
247 263 /* u+v = s*(1+...) */
248 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ 264 u = s_h * t_h;
249 if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0) 265 v = s_l * t_h + t_l * s;
250 s = -one;/* (-ve)**(odd int) */ 266 /* 2/(3log2)*(s+...) */
267 p_h = u + v;
268 SET_LOW_WORD(p_h, 0);
269 p_l = v - (p_h - u);
270 z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
271 z_l = cp_l * p_h + p_l * cp + dp_l[k];
272 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
273 t = (double) n;
274 t1 = (((z_h + z_l) + dp_h[k]) + t);
275 SET_LOW_WORD(t1, 0);
276 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
277 }
278
279 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
280 if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
281 s = -one; /* (-ve)**(odd int) */
251 282
252 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ 283 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
253 y1 = y; 284 y1 = y;
254 SET_LOW_WORD(y1,0); 285 SET_LOW_WORD(y1, 0);
255 p_l = (y-y1)*t1+y*t2; 286 p_l = (y - y1) * t1 + y * t2;
256 p_h = y1*t1; 287 p_h = y1 * t1;
257 z = p_l+p_h; 288 z = p_l + p_h;
258 EXTRACT_WORDS(j,i,z); 289 EXTRACT_WORDS(j, i, z);
259 if (j>=0x40900000) { /* z >= 1024 */ 290 if (j >= 0x40900000) { /* z >= 1024 */
260 if(((j-0x40900000)|i)!=0) /* if z > 1024 */ 291 if (((j - 0x40900000) | i) != 0) /* if z > 1024 */
261 return s*huge*huge; /* overflow */ 292 return s * huge * huge; /* overflow */
262 else { 293 else {
263 if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ 294 if (p_l + ovt > z - p_h)
264 } 295 return s * huge * huge; /* overflow */
265 } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */ 296 }
266 if(((j-0xc090cc00)|i)!=0) /* z < -1075 */ 297 } else if ((j & 0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */
267 return s*tiny*tiny; /* underflow */ 298 if (((j - 0xc090cc00) | i) != 0) /* z < -1075 */
268 else { 299 return s * tiny * tiny; /* underflow */
269 if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ 300 else {
270 } 301 if (p_l <= z - p_h)
271 } 302 return s * tiny * tiny; /* underflow */
303 }
304 }
272 /* 305 /*
273 * compute 2**(p_h+p_l) 306 * compute 2**(p_h+p_l)
274 */ 307 */
275 i = j&0x7fffffff; 308 i = j & 0x7fffffff;
276 k = (i>>20)-0x3ff; 309 k = (i >> 20) - 0x3ff;
277 n = 0; 310 n = 0;
278 if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ 311 if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
279 n = j+(0x00100000>>(k+1)); 312 n = j + (0x00100000 >> (k + 1));
280 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */ 313 k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */
281 t = zero; 314 t = zero;
282 SET_HIGH_WORD(t,n&~(0x000fffff>>k)); 315 SET_HIGH_WORD(t, n & ~(0x000fffff >> k));
283 n = ((n&0x000fffff)|0x00100000)>>(20-k); 316 n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
284 if(j<0) n = -n; 317 if (j < 0)
285 p_h -= t; 318 n = -n;
286 } 319 p_h -= t;
287 t = p_l+p_h; 320 }
288 SET_LOW_WORD(t,0); 321 t = p_l + p_h;
289 u = t*lg2_h; 322 SET_LOW_WORD(t, 0);
290 v = (p_l-(t-p_h))*lg2+t*lg2_l; 323 u = t * lg2_h;
291 z = u+v; 324 v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
292 w = v-(z-u); 325 z = u + v;
293 t = z*z; 326 w = v - (z - u);
294 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 327 t = z * z;
295 r = (z*t1)/(t1-two)-(w+z*w); 328 t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
296 z = one-(r-z); 329 r = (z * t1) / (t1 - two) - (w + z * w);
297 GET_HIGH_WORD(j,z); 330 z = one - (r - z);
298 j += (n<<20); 331 GET_HIGH_WORD(j, z);
299 if((j>>20)<=0) z = SDL_NAME(scalbn)(z,n); /* subnormal output */ 332 j += (n << 20);
300 else SET_HIGH_WORD(z,j); 333 if ((j >> 20) <= 0)
301 return s*z; 334 z = SDL_NAME(scalbn) (z, n); /* subnormal output */
335 else
336 SET_HIGH_WORD(z, j);
337 return s * z;
302 } 338 }
339
340 /* vi: set ts=4 sw=4 expandtab: */