- /*
- * We no longer need to avoid falling into the multi-precision
- * calculations due to compiler bugs breaking Dekker's theorem.
- * Keep avoiding this as an optimization. See log2.c for more
- * details (some details are here only because the optimization
- * is not yet available in double precision).
- *
- * Another compiler bug turned up. With gcc on i386,
- * (ivln2lo + ivln2hi) would be evaluated in float precision
- * despite runtime evaluations using double precision. So we
- * must cast one of its terms to float_t. This makes the whole
- * expression have type float_t, so return is forced to waste
- * time clobbering its extra precision.
- */
-// FIXME
-// if (sizeof(float_t) > sizeof(float))
-// return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y;
+ /* reduce x into [sqrt(2)/2, sqrt(2)] */
+ ix += 0x3f800000 - 0x3f3504f3;
+ k += (int)(ix>>23) - 0x7f;
+ ix = (ix&0x007fffff) + 0x3f3504f3;
+ u.i = ix;
+ x = u.f;
+
+ f = x - 1.0f;
+ s = f/(2.0f + f);
+ z = s*s;
+ w = z*z;
+ t1= w*(Lg2+w*Lg4);
+ t2= z*(Lg1+w*Lg3);
+ R = t2 + t1;
+ hfsq = 0.5f*f*f;