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

Side by Side Diff: fusl/src/math/log10.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_log10.c */ 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log10.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 10 logarithm of x. See log.c for most comments. 13 * Return the base 10 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 * log10(x) = (f - f*f/2 + r)/log(10) + k*log10(2) 17 * log10(x) = (f - f*f/2 + r)/log(10) + k*log10(2)
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 ivln10hi =
24 ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */ 24 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */
25 ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ 25 ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
26 log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ 26 log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
27 log10_2lo = 3.69423907715893078616e-13, /* 0x3D59FEF3, 0x11F12B36 */ 27 log10_2lo = 3.69423907715893078616e-13, /* 0x3D59FEF3, 0x11F12B36 */
28 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 28 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
29 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 29 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
30 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 30 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
31 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 31 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
32 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 32 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
33 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 33 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
34 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 34 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
35 35
36 double log10(double x) 36 double log10(double x) {
37 { 37 union {
38 » union {double f; uint64_t i;} u = {x}; 38 double f;
39 » double_t hfsq,f,s,z,R,w,t1,t2,dk,y,hi,lo,val_hi,val_lo; 39 uint64_t i;
40 » uint32_t hx; 40 } u = {x};
41 » int k; 41 double_t hfsq, f, s, z, R, w, t1, t2, dk, y, hi, lo, val_hi, val_lo;
42 uint32_t hx;
43 int k;
42 44
43 » hx = u.i>>32; 45 hx = u.i >> 32;
44 » k = 0; 46 k = 0;
45 » if (hx < 0x00100000 || hx>>31) { 47 if (hx < 0x00100000 || hx >> 31) {
46 » » if (u.i<<1 == 0) 48 if (u.i << 1 == 0)
47 » » » return -1/(x*x); /* log(+-0)=-inf */ 49 return -1 / (x * x); /* log(+-0)=-inf */
48 » » if (hx>>31) 50 if (hx >> 31)
49 » » » return (x-x)/0.0; /* log(-#) = NaN */ 51 return (x - x) / 0.0; /* log(-#) = NaN */
50 » » /* subnormal number, scale x up */ 52 /* subnormal number, scale x up */
51 » » k -= 54; 53 k -= 54;
52 » » x *= 0x1p54; 54 x *= 0x1p54;
53 » » u.f = x; 55 u.f = x;
54 » » hx = u.i>>32; 56 hx = u.i >> 32;
55 » } else if (hx >= 0x7ff00000) { 57 } else if (hx >= 0x7ff00000) {
56 » » return x; 58 return x;
57 » } else if (hx == 0x3ff00000 && u.i<<32 == 0) 59 } else if (hx == 0x3ff00000 && u.i << 32 == 0)
58 » » return 0; 60 return 0;
59 61
60 » /* reduce x into [sqrt(2)/2, sqrt(2)] */ 62 /* reduce x into [sqrt(2)/2, sqrt(2)] */
61 » hx += 0x3ff00000 - 0x3fe6a09e; 63 hx += 0x3ff00000 - 0x3fe6a09e;
62 » k += (int)(hx>>20) - 0x3ff; 64 k += (int)(hx >> 20) - 0x3ff;
63 » hx = (hx&0x000fffff) + 0x3fe6a09e; 65 hx = (hx & 0x000fffff) + 0x3fe6a09e;
64 » u.i = (uint64_t)hx<<32 | (u.i&0xffffffff); 66 u.i = (uint64_t)hx << 32 | (u.i & 0xffffffff);
65 » x = u.f; 67 x = u.f;
66 68
67 » f = x - 1.0; 69 f = x - 1.0;
68 » hfsq = 0.5*f*f; 70 hfsq = 0.5 * f * f;
69 » s = f/(2.0+f); 71 s = f / (2.0 + f);
70 » z = s*s; 72 z = s * s;
71 » w = z*z; 73 w = z * z;
72 » t1 = w*(Lg2+w*(Lg4+w*Lg6)); 74 t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
73 » t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 75 t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
74 » R = t2 + t1; 76 R = t2 + t1;
75 77
76 » /* See log2.c for details. */ 78 /* See log2.c for details. */
77 » /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */ 79 /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
78 » hi = f - hfsq; 80 hi = f - hfsq;
79 » u.f = hi; 81 u.f = hi;
80 » u.i &= (uint64_t)-1<<32; 82 u.i &= (uint64_t)-1 << 32;
81 » hi = u.f; 83 hi = u.f;
82 » lo = f - hi - hfsq + s*(hfsq+R); 84 lo = f - hi - hfsq + s * (hfsq + R);
83 85
84 » /* val_hi+val_lo ~ log10(1+f) + k*log10(2) */ 86 /* val_hi+val_lo ~ log10(1+f) + k*log10(2) */
85 » val_hi = hi*ivln10hi; 87 val_hi = hi * ivln10hi;
86 » dk = k; 88 dk = k;
87 » y = dk*log10_2hi; 89 y = dk * log10_2hi;
88 » val_lo = dk*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi; 90 val_lo = dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi;
89 91
90 » /* 92 /*
91 » * Extra precision in for adding y is not strictly needed 93 * Extra precision in for adding y is not strictly needed
92 » * since there is no very large cancellation near x = sqrt(2) or 94 * since there is no very large cancellation near x = sqrt(2) or
93 » * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs 95 * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
94 » * with some parallelism and it reduces the error for many args. 96 * with some parallelism and it reduces the error for many args.
95 » */ 97 */
96 » w = y + val_hi; 98 w = y + val_hi;
97 » val_lo += (y - w) + val_hi; 99 val_lo += (y - w) + val_hi;
98 » val_hi = w; 100 val_hi = w;
99 101
100 » return val_lo + val_hi; 102 return val_lo + val_hi;
101 } 103 }
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698