Branch data Line data Source code
1 : :
2 : : /* @(#)e_log10.c 1.3 95/01/18 */
3 : : /*
4 : : * ====================================================
5 : : * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 : : *
7 : : * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 : : * Permission to use, copy, modify, and distribute this
9 : : * software is freely granted, provided that this notice
10 : : * is preserved.
11 : : * ====================================================
12 : : */
13 : :
14 : : #include "cdefs-compat.h"
15 : : //__FBSDID("$FreeBSD: src/lib/msun/src/e_log2.c,v 1.4 2011/10/15 05:23:28 das Exp $");
16 : :
17 : : /*
18 : : * Return the base 2 logarithm of x. See e_log.c and k_log.h for most
19 : : * comments.
20 : : *
21 : : * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
22 : : * then does the combining and scaling steps
23 : : * log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
24 : : * in not-quite-routine extra precision.
25 : : */
26 : :
27 : : #include <float.h>
28 : : #include <openlibm_math.h>
29 : :
30 : : #include "math_private.h"
31 : : #include "k_log.h"
32 : :
33 : : static const double
34 : : two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
35 : : ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
36 : : ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
37 : :
38 : : static const double zero = 0.0;
39 : :
40 : : OLM_DLLEXPORT double
41 : 12 : __ieee754_log2(double x)
42 : : {
43 : : double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
44 : : int32_t i,k,hx;
45 : : u_int32_t lx;
46 : :
47 : 12 : EXTRACT_WORDS(hx,lx,x);
48 : :
49 : 12 : k=0;
50 [ + + ]: 12 : if (hx < 0x00100000) { /* x < 2**-1022 */
51 [ + + ]: 4 : if (((hx&0x7fffffff)|lx)==0)
52 : 2 : return -two54/zero; /* log(+-0)=-inf */
53 [ + - ]: 2 : if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
54 : 0 : k -= 54; x *= two54; /* subnormal number, scale up x */
55 : 0 : GET_HIGH_WORD(hx,x);
56 : : }
57 [ + + ]: 8 : if (hx >= 0x7ff00000) return x+x;
58 [ + + + - ]: 7 : if (hx == 0x3ff00000 && lx == 0)
59 : 2 : return zero; /* log(1) = +0 */
60 : 5 : k += (hx>>20)-1023;
61 : 5 : hx &= 0x000fffff;
62 : 5 : i = (hx+0x95f64)&0x100000;
63 : 5 : SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */
64 : 5 : k += (i>>20);
65 : 5 : y = (double)k;
66 : 5 : f = x - 1.0;
67 : 5 : hfsq = 0.5*f*f;
68 : 5 : r = k_log1p(f);
69 : :
70 : : /*
71 : : * f-hfsq must (for args near 1) be evaluated in extra precision
72 : : * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
73 : : * This is fairly efficient since f-hfsq only depends on f, so can
74 : : * be evaluated in parallel with R. Not combining hfsq with R also
75 : : * keeps R small (though not as small as a true `lo' term would be),
76 : : * so that extra precision is not needed for terms involving R.
77 : : *
78 : : * Compiler bugs involving extra precision used to break Dekker's
79 : : * theorem for spitting f-hfsq as hi+lo, unless double_t was used
80 : : * or the multi-precision calculations were avoided when double_t
81 : : * has extra precision. These problems are now automatically
82 : : * avoided as a side effect of the optimization of combining the
83 : : * Dekker splitting step with the clear-low-bits step.
84 : : *
85 : : * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
86 : : * precision to avoid a very large cancellation when x is very near
87 : : * these values. Unlike the above cancellations, this problem is
88 : : * specific to base 2. It is strange that adding +-1 is so much
89 : : * harder than adding +-ln2 or +-log10_2.
90 : : *
91 : : * This uses Dekker's theorem to normalize y+val_hi, so the
92 : : * compiler bugs are back in some configurations, sigh. And I
93 : : * don't want to used double_t to avoid them, since that gives a
94 : : * pessimization and the support for avoiding the pessimization
95 : : * is not yet available.
96 : : *
97 : : * The multi-precision calculations for the multiplications are
98 : : * routine.
99 : : */
100 : 5 : hi = f - hfsq;
101 : 5 : SET_LOW_WORD(hi,0);
102 : 5 : lo = (f - hi) - hfsq + r;
103 : 5 : val_hi = hi*ivln2hi;
104 : 5 : val_lo = (lo+hi)*ivln2lo + lo*ivln2hi;
105 : :
106 : : /* spadd(val_hi, val_lo, y), except for not using double_t: */
107 : 5 : w = y + val_hi;
108 : 5 : val_lo += (y - w) + val_hi;
109 : 5 : val_hi = w;
110 : :
111 : 5 : return val_lo + val_hi;
112 : : }
113 : :
114 : : #if (LDBL_MANT_DIG == 53)
115 : : openlibm_weak_reference(log2, log2l);
116 : : #endif
|
