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

Side by Side Diff: fusl/src/math/log1p.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/s_log1p.c */ 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.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 SunPro, a Sun Microsystems, Inc. business. 6 * Developed at SunPro, 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 * ====================================================
(...skipping 37 matching lines...) Expand 10 before | Expand all | Expand 10 after
48 * 48 *
49 * u = 1+x; 49 * u = 1+x;
50 * if(u==1.0) return x ; else 50 * if(u==1.0) return x ; else
51 * return log(u)*(x/(u-1.0)); 51 * return log(u)*(x/(u-1.0));
52 * 52 *
53 * See HP-15C Advanced Functions Handbook, p.193. 53 * See HP-15C Advanced Functions Handbook, p.193.
54 */ 54 */
55 55
56 #include "libm.h" 56 #include "libm.h"
57 57
58 static const double 58 static const double ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
59 ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ 59 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
60 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ 60 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
61 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 61 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
62 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 62 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
63 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 63 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
64 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 64 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
65 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 65 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
66 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 66 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
67 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
68 67
69 double log1p(double x) 68 double log1p(double x) {
70 { 69 union {
71 » union {double f; uint64_t i;} u = {x}; 70 double f;
72 » double_t hfsq,f,c,s,z,R,w,t1,t2,dk; 71 uint64_t i;
73 » uint32_t hx,hu; 72 } u = {x};
74 » int k; 73 double_t hfsq, f, c, s, z, R, w, t1, t2, dk;
74 uint32_t hx, hu;
75 int k;
75 76
76 » hx = u.i>>32; 77 hx = u.i >> 32;
77 » k = 1; 78 k = 1;
78 » if (hx < 0x3fda827a || hx>>31) { /* 1+x < sqrt(2)+ */ 79 if (hx < 0x3fda827a || hx >> 31) { /* 1+x < sqrt(2)+ */
79 » » if (hx >= 0xbff00000) { /* x <= -1.0 */ 80 if (hx >= 0xbff00000) { /* x <= -1.0 */
80 » » » if (x == -1) 81 if (x == -1)
81 » » » » return x/0.0; /* log1p(-1) = -inf */ 82 return x / 0.0; /* log1p(-1) = -inf */
82 » » » return (x-x)/0.0; /* log1p(x<-1) = NaN */ 83 return (x - x) / 0.0; /* log1p(x<-1) = NaN */
83 » » } 84 }
84 » » if (hx<<1 < 0x3ca00000<<1) { /* |x| < 2**-53 */ 85 if (hx << 1 < 0x3ca00000 << 1) { /* |x| < 2**-53 */
85 » » » /* underflow if subnormal */ 86 /* underflow if subnormal */
86 » » » if ((hx&0x7ff00000) == 0) 87 if ((hx & 0x7ff00000) == 0)
87 » » » » FORCE_EVAL((float)x); 88 FORCE_EVAL((float)x);
88 » » » return x; 89 return x;
89 » » } 90 }
90 » » if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ 91 if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
91 » » » k = 0; 92 k = 0;
92 » » » c = 0; 93 c = 0;
93 » » » f = x; 94 f = x;
94 » » } 95 }
95 » } else if (hx >= 0x7ff00000) 96 } else if (hx >= 0x7ff00000)
96 » » return x; 97 return x;
97 » if (k) { 98 if (k) {
98 » » u.f = 1 + x; 99 u.f = 1 + x;
99 » » hu = u.i>>32; 100 hu = u.i >> 32;
100 » » hu += 0x3ff00000 - 0x3fe6a09e; 101 hu += 0x3ff00000 - 0x3fe6a09e;
101 » » k = (int)(hu>>20) - 0x3ff; 102 k = (int)(hu >> 20) - 0x3ff;
102 » » /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ 103 /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
103 » » if (k < 54) { 104 if (k < 54) {
104 » » » c = k >= 2 ? 1-(u.f-x) : x-(u.f-1); 105 c = k >= 2 ? 1 - (u.f - x) : x - (u.f - 1);
105 » » » c /= u.f; 106 c /= u.f;
106 » » } else 107 } else
107 » » » c = 0; 108 c = 0;
108 » » /* reduce u into [sqrt(2)/2, sqrt(2)] */ 109 /* reduce u into [sqrt(2)/2, sqrt(2)] */
109 » » hu = (hu&0x000fffff) + 0x3fe6a09e; 110 hu = (hu & 0x000fffff) + 0x3fe6a09e;
110 » » u.i = (uint64_t)hu<<32 | (u.i&0xffffffff); 111 u.i = (uint64_t)hu << 32 | (u.i & 0xffffffff);
111 » » f = u.f - 1; 112 f = u.f - 1;
112 » } 113 }
113 » hfsq = 0.5*f*f; 114 hfsq = 0.5 * f * f;
114 » s = f/(2.0+f); 115 s = f / (2.0 + f);
115 » z = s*s; 116 z = s * s;
116 » w = z*z; 117 w = z * z;
117 » t1 = w*(Lg2+w*(Lg4+w*Lg6)); 118 t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
118 » t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 119 t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
119 » R = t2 + t1; 120 R = t2 + t1;
120 » dk = k; 121 dk = k;
121 » return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi; 122 return s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi;
122 } 123 }
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698