Chromium Code Reviews
chromiumcodereview-hr@appspot.gserviceaccount.com (chromiumcodereview-hr) | Please choose your nickname with Settings | Help | Chromium Project | Gerrit Changes | Sign out
(132)

Side by Side Diff: fusl/src/math/log2.c

Issue 1714623002: [fusl] clang-format fusl (Closed) Base URL: git@github.com:domokit/mojo.git@master
Patch Set: headers too Created 4 years, 10 months ago
Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.
Jump to:
View unified diff | Download patch
OLDNEW
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log2.c */ 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log2.c */
2 /* 2 /*
3 * ==================================================== 3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 * 5 *
6 * Developed at SunSoft, a Sun Microsystems, Inc. business. 6 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this 7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice 8 * software is freely granted, provided that this notice
9 * is preserved. 9 * is preserved.
10 * ==================================================== 10 * ====================================================
11 */ 11 */
12 /* 12 /*
13 * Return the base 2 logarithm of x. See log.c for most comments. 13 * Return the base 2 logarithm of x. See log.c for most comments.
14 * 14 *
15 * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2 15 * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2
16 * as in log.c, then combine and scale in extra precision: 16 * as in log.c, then combine and scale in extra precision:
17 * log2(x) = (f - f*f/2 + r)/log(2) + k 17 * log2(x) = (f - f*f/2 + r)/log(2) + k
18 */ 18 */
19 19
20 #include <math.h> 20 #include <math.h>
21 #include <stdint.h> 21 #include <stdint.h>
22 22
23 static const double 23 static const double ivln2hi =
24 ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */ 24 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
25 ivln2lo = 1.67517131648865118353e-10, /* 0x3de705fc, 0x2eefa200 */ 25 ivln2lo = 1.67517131648865118353e-10, /* 0x3de705fc, 0x2eefa200 */
26 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 26 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
27 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 27 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
28 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 28 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
29 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 29 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
30 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 30 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
31 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 31 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
32 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 32 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
33 33
34 double log2(double x) 34 double log2(double x) {
35 { 35 union {
36 » union {double f; uint64_t i;} u = {x}; 36 double f;
37 » double_t hfsq,f,s,z,R,w,t1,t2,y,hi,lo,val_hi,val_lo; 37 uint64_t i;
38 » uint32_t hx; 38 } u = {x};
39 » int k; 39 double_t hfsq, f, s, z, R, w, t1, t2, y, hi, lo, val_hi, val_lo;
40 uint32_t hx;
41 int k;
40 42
41 » hx = u.i>>32; 43 hx = u.i >> 32;
42 » k = 0; 44 k = 0;
43 » if (hx < 0x00100000 || hx>>31) { 45 if (hx < 0x00100000 || hx >> 31) {
44 » » if (u.i<<1 == 0) 46 if (u.i << 1 == 0)
45 » » » return -1/(x*x); /* log(+-0)=-inf */ 47 return -1 / (x * x); /* log(+-0)=-inf */
46 » » if (hx>>31) 48 if (hx >> 31)
47 » » » return (x-x)/0.0; /* log(-#) = NaN */ 49 return (x - x) / 0.0; /* log(-#) = NaN */
48 » » /* subnormal number, scale x up */ 50 /* subnormal number, scale x up */
49 » » k -= 54; 51 k -= 54;
50 » » x *= 0x1p54; 52 x *= 0x1p54;
51 » » u.f = x; 53 u.f = x;
52 » » hx = u.i>>32; 54 hx = u.i >> 32;
53 » } else if (hx >= 0x7ff00000) { 55 } else if (hx >= 0x7ff00000) {
54 » » return x; 56 return x;
55 » } else if (hx == 0x3ff00000 && u.i<<32 == 0) 57 } else if (hx == 0x3ff00000 && u.i << 32 == 0)
56 » » return 0; 58 return 0;
57 59
58 » /* reduce x into [sqrt(2)/2, sqrt(2)] */ 60 /* reduce x into [sqrt(2)/2, sqrt(2)] */
59 » hx += 0x3ff00000 - 0x3fe6a09e; 61 hx += 0x3ff00000 - 0x3fe6a09e;
60 » k += (int)(hx>>20) - 0x3ff; 62 k += (int)(hx >> 20) - 0x3ff;
61 » hx = (hx&0x000fffff) + 0x3fe6a09e; 63 hx = (hx & 0x000fffff) + 0x3fe6a09e;
62 » u.i = (uint64_t)hx<<32 | (u.i&0xffffffff); 64 u.i = (uint64_t)hx << 32 | (u.i & 0xffffffff);
63 » x = u.f; 65 x = u.f;
64 66
65 » f = x - 1.0; 67 f = x - 1.0;
66 » hfsq = 0.5*f*f; 68 hfsq = 0.5 * f * f;
67 » s = f/(2.0+f); 69 s = f / (2.0 + f);
68 » z = s*s; 70 z = s * s;
69 » w = z*z; 71 w = z * z;
70 » t1 = w*(Lg2+w*(Lg4+w*Lg6)); 72 t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
71 » t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 73 t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
72 » R = t2 + t1; 74 R = t2 + t1;
73 75
74 » /* 76 /*
75 » * f-hfsq must (for args near 1) be evaluated in extra precision 77 * f-hfsq must (for args near 1) be evaluated in extra precision
76 » * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2). 78 * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
77 » * This is fairly efficient since f-hfsq only depends on f, so can 79 * This is fairly efficient since f-hfsq only depends on f, so can
78 » * be evaluated in parallel with R. Not combining hfsq with R also 80 * be evaluated in parallel with R. Not combining hfsq with R also
79 » * keeps R small (though not as small as a true `lo' term would be), 81 * keeps R small (though not as small as a true `lo' term would be),
80 » * so that extra precision is not needed for terms involving R. 82 * so that extra precision is not needed for terms involving R.
81 » * 83 *
82 » * Compiler bugs involving extra precision used to break Dekker's 84 * Compiler bugs involving extra precision used to break Dekker's
83 » * theorem for spitting f-hfsq as hi+lo, unless double_t was used 85 * theorem for spitting f-hfsq as hi+lo, unless double_t was used
84 » * or the multi-precision calculations were avoided when double_t 86 * or the multi-precision calculations were avoided when double_t
85 » * has extra precision. These problems are now automatically 87 * has extra precision. These problems are now automatically
86 » * avoided as a side effect of the optimization of combining the 88 * avoided as a side effect of the optimization of combining the
87 » * Dekker splitting step with the clear-low-bits step. 89 * Dekker splitting step with the clear-low-bits step.
88 » * 90 *
89 » * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra 91 * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
90 » * precision to avoid a very large cancellation when x is very near 92 * precision to avoid a very large cancellation when x is very near
91 » * these values. Unlike the above cancellations, this problem is 93 * these values. Unlike the above cancellations, this problem is
92 » * specific to base 2. It is strange that adding +-1 is so much 94 * specific to base 2. It is strange that adding +-1 is so much
93 » * harder than adding +-ln2 or +-log10_2. 95 * harder than adding +-ln2 or +-log10_2.
94 » * 96 *
95 » * This uses Dekker's theorem to normalize y+val_hi, so the 97 * This uses Dekker's theorem to normalize y+val_hi, so the
96 » * compiler bugs are back in some configurations, sigh. And I 98 * compiler bugs are back in some configurations, sigh. And I
97 » * don't want to used double_t to avoid them, since that gives a 99 * don't want to used double_t to avoid them, since that gives a
98 » * pessimization and the support for avoiding the pessimization 100 * pessimization and the support for avoiding the pessimization
99 » * is not yet available. 101 * is not yet available.
100 » * 102 *
101 » * The multi-precision calculations for the multiplications are 103 * The multi-precision calculations for the multiplications are
102 » * routine. 104 * routine.
103 » */ 105 */
104 106
105 » /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */ 107 /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
106 » hi = f - hfsq; 108 hi = f - hfsq;
107 » u.f = hi; 109 u.f = hi;
108 » u.i &= (uint64_t)-1<<32; 110 u.i &= (uint64_t)-1 << 32;
109 » hi = u.f; 111 hi = u.f;
110 » lo = f - hi - hfsq + s*(hfsq+R); 112 lo = f - hi - hfsq + s * (hfsq + R);
111 113
112 » val_hi = hi*ivln2hi; 114 val_hi = hi * ivln2hi;
113 » val_lo = (lo+hi)*ivln2lo + lo*ivln2hi; 115 val_lo = (lo + hi) * ivln2lo + lo * ivln2hi;
114 116
115 » /* spadd(val_hi, val_lo, y), except for not using double_t: */ 117 /* spadd(val_hi, val_lo, y), except for not using double_t: */
116 » y = k; 118 y = k;
117 » w = y + val_hi; 119 w = y + val_hi;
118 » val_lo += (y - w) + val_hi; 120 val_lo += (y - w) + val_hi;
119 » val_hi = w; 121 val_hi = w;
120 122
121 » return val_lo + val_hi; 123 return val_lo + val_hi;
122 } 124 }
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698