1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log2f.c */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
10 * ====================================================
13 * See comments in log2.c.
20 two25 = 3.3554432000e+07, /* 0x4c000000 */
21 ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */
22 ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
24 static const float zero = 0.0;
28 float f,hfsq,hi,lo,r,y;
31 GET_FLOAT_WORD(hx, x);
34 if (hx < 0x00800000) { /* x < 2**-126 */
35 if ((hx&0x7fffffff) == 0)
36 return -two25/zero; /* log(+-0)=-inf */
38 return (x-x)/zero; /* log(-#) = NaN */
39 /* subnormal number, scale up x */
42 GET_FLOAT_WORD(hx, x);
47 return zero; /* log(1) = +0 */
50 i = (hx+(0x4afb0d))&0x800000;
51 SET_FLOAT_WORD(x, hx|(i^0x3f800000)); /* normalize x or x/2 */
59 * We no longer need to avoid falling into the multi-precision
60 * calculations due to compiler bugs breaking Dekker's theorem.
61 * Keep avoiding this as an optimization. See log2.c for more
62 * details (some details are here only because the optimization
63 * is not yet available in double precision).
65 * Another compiler bug turned up. With gcc on i386,
66 * (ivln2lo + ivln2hi) would be evaluated in float precision
67 * despite runtime evaluations using double precision. So we
68 * must cast one of its terms to float_t. This makes the whole
69 * expression have type float_t, so return is forced to waste
70 * time clobbering its extra precision.
73 // if (sizeof(float_t) > sizeof(float))
74 // return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y;
77 GET_FLOAT_WORD(hx,hi);
78 SET_FLOAT_WORD(hi,hx&0xfffff000);
79 lo = (f - hi) - hfsq + r;
80 return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + y;