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

Side by Side Diff: fusl/src/math/log1pl.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: OpenBSD /usr/src/lib/libm/src/ld80/s_log1pl.c */ 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_log1pl.c */
2 /* 2 /*
3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
4 * 4 *
5 * Permission to use, copy, modify, and distribute this software for any 5 * Permission to use, copy, modify, and distribute this software for any
6 * purpose with or without fee is hereby granted, provided that the above 6 * purpose with or without fee is hereby granted, provided that the above
7 * copyright notice and this permission notice appear in all copies. 7 * copyright notice and this permission notice appear in all copies.
8 * 8 *
9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
(...skipping 33 matching lines...) Expand 10 before | Expand all | Expand 10 after
44 * ACCURACY: 44 * ACCURACY:
45 * 45 *
46 * Relative error: 46 * Relative error:
47 * arithmetic domain # trials peak rms 47 * arithmetic domain # trials peak rms
48 * IEEE -1.0, 9.0 100000 8.2e-20 2.5e-20 48 * IEEE -1.0, 9.0 100000 8.2e-20 2.5e-20
49 */ 49 */
50 50
51 #include "libm.h" 51 #include "libm.h"
52 52
53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
54 long double log1pl(long double x) 54 long double log1pl(long double x) {
55 { 55 return log1p(x);
56 » return log1p(x);
57 } 56 }
58 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 57 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
59 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x) 58 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
60 * 1/sqrt(2) <= x < sqrt(2) 59 * 1/sqrt(2) <= x < sqrt(2)
61 * Theoretical peak relative error = 2.32e-20 60 * Theoretical peak relative error = 2.32e-20
62 */ 61 */
63 static const long double P[] = { 62 static const long double P[] = {
64 4.5270000862445199635215E-5L, 63 4.5270000862445199635215E-5L, 4.9854102823193375972212E-1L,
65 4.9854102823193375972212E-1L, 64 6.5787325942061044846969E0L, 2.9911919328553073277375E1L,
66 6.5787325942061044846969E0L, 65 6.0949667980987787057556E1L, 5.7112963590585538103336E1L,
67 2.9911919328553073277375E1L, 66 2.0039553499201281259648E1L,
68 6.0949667980987787057556E1L,
69 5.7112963590585538103336E1L,
70 2.0039553499201281259648E1L,
71 }; 67 };
72 static const long double Q[] = { 68 static const long double Q[] = {
73 /* 1.0000000000000000000000E0,*/ 69 /* 1.0000000000000000000000E0,*/
74 1.5062909083469192043167E1L, 70 1.5062909083469192043167E1L, 8.3047565967967209469434E1L,
75 8.3047565967967209469434E1L, 71 2.2176239823732856465394E2L, 3.0909872225312059774938E2L,
76 2.2176239823732856465394E2L, 72 2.1642788614495947685003E2L, 6.0118660497603843919306E1L,
77 3.0909872225312059774938E2L,
78 2.1642788614495947685003E2L,
79 6.0118660497603843919306E1L,
80 }; 73 };
81 74
82 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), 75 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
83 * where z = 2(x-1)/(x+1) 76 * where z = 2(x-1)/(x+1)
84 * 1/sqrt(2) <= x < sqrt(2) 77 * 1/sqrt(2) <= x < sqrt(2)
85 * Theoretical peak relative error = 6.16e-22 78 * Theoretical peak relative error = 6.16e-22
86 */ 79 */
87 static const long double R[4] = { 80 static const long double R[4] = {
88 1.9757429581415468984296E-3L, 81 1.9757429581415468984296E-3L, -7.1990767473014147232598E-1L,
89 -7.1990767473014147232598E-1L, 82 1.0777257190312272158094E1L, -3.5717684488096787370998E1L,
90 1.0777257190312272158094E1L,
91 -3.5717684488096787370998E1L,
92 }; 83 };
93 static const long double S[4] = { 84 static const long double S[4] = {
94 /* 1.00000000000000000000E0L,*/ 85 /* 1.00000000000000000000E0L,*/
95 -2.6201045551331104417768E1L, 86 -2.6201045551331104417768E1L, 1.9361891836232102174846E2L,
96 1.9361891836232102174846E2L, 87 -4.2861221385716144629696E2L,
97 -4.2861221385716144629696E2L,
98 }; 88 };
99 static const long double C1 = 6.9314575195312500000000E-1L; 89 static const long double C1 = 6.9314575195312500000000E-1L;
100 static const long double C2 = 1.4286068203094172321215E-6L; 90 static const long double C2 = 1.4286068203094172321215E-6L;
101 91
102 #define SQRTH 0.70710678118654752440L 92 #define SQRTH 0.70710678118654752440L
103 93
104 long double log1pl(long double xm1) 94 long double log1pl(long double xm1) {
105 { 95 long double x, y, z;
106 » long double x, y, z; 96 int e;
107 » int e;
108 97
109 » if (isnan(xm1)) 98 if (isnan(xm1))
110 » » return xm1; 99 return xm1;
111 » if (xm1 == INFINITY) 100 if (xm1 == INFINITY)
112 » » return xm1; 101 return xm1;
113 » if (xm1 == 0.0) 102 if (xm1 == 0.0)
114 » » return xm1; 103 return xm1;
115 104
116 » x = xm1 + 1.0; 105 x = xm1 + 1.0;
117 106
118 » /* Test for domain errors. */ 107 /* Test for domain errors. */
119 » if (x <= 0.0) { 108 if (x <= 0.0) {
120 » » if (x == 0.0) 109 if (x == 0.0)
121 » » » return -1/(x*x); /* -inf with divbyzero */ 110 return -1 / (x * x); /* -inf with divbyzero */
122 » » return 0/0.0f; /* nan with invalid */ 111 return 0 / 0.0f; /* nan with invalid */
123 » } 112 }
124 113
125 » /* Separate mantissa from exponent. 114 /* Separate mantissa from exponent.
126 » Use frexp so that denormal numbers will be handled properly. */ 115 Use frexp so that denormal numbers will be handled properly. */
127 » x = frexpl(x, &e); 116 x = frexpl(x, &e);
128 117
129 » /* logarithm using log(x) = z + z^3 P(z)/Q(z), 118 /* logarithm using log(x) = z + z^3 P(z)/Q(z),
130 » where z = 2(x-1)/x+1) */ 119 where z = 2(x-1)/x+1) */
131 » if (e > 2 || e < -2) { 120 if (e > 2 || e < -2) {
132 » » if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ 121 if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
133 » » » e -= 1; 122 e -= 1;
134 » » » z = x - 0.5; 123 z = x - 0.5;
135 » » » y = 0.5 * z + 0.5; 124 y = 0.5 * z + 0.5;
136 » » } else { /* 2 (x-1)/(x+1) */ 125 } else { /* 2 (x-1)/(x+1) */
137 » » » z = x - 0.5; 126 z = x - 0.5;
138 » » » z -= 0.5; 127 z -= 0.5;
139 » » » y = 0.5 * x + 0.5; 128 y = 0.5 * x + 0.5;
140 » » } 129 }
141 » » x = z / y; 130 x = z / y;
142 » » z = x*x; 131 z = x * x;
143 » » z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); 132 z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3));
144 » » z = z + e * C2; 133 z = z + e * C2;
145 » » z = z + x; 134 z = z + x;
146 » » z = z + e * C1; 135 z = z + e * C1;
147 » » return z; 136 return z;
148 » } 137 }
149 138
150 » /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ 139 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
151 » if (x < SQRTH) { 140 if (x < SQRTH) {
152 » » e -= 1; 141 e -= 1;
153 » » if (e != 0) 142 if (e != 0)
154 » » » x = 2.0 * x - 1.0; 143 x = 2.0 * x - 1.0;
155 » » else 144 else
156 » » » x = xm1; 145 x = xm1;
157 » } else { 146 } else {
158 » » if (e != 0) 147 if (e != 0)
159 » » » x = x - 1.0; 148 x = x - 1.0;
160 » » else 149 else
161 » » » x = xm1; 150 x = xm1;
162 » } 151 }
163 » z = x*x; 152 z = x * x;
164 » y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6)); 153 y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
165 » y = y + e * C2; 154 y = y + e * C2;
166 » z = y - 0.5 * z; 155 z = y - 0.5 * z;
167 » z = z + x; 156 z = z + x;
168 » z = z + e * C1; 157 z = z + e * C1;
169 » return z; 158 return z;
170 } 159 }
171 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 160 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
172 // TODO: broken implementation to make things compile 161 // TODO: broken implementation to make things compile
173 long double log1pl(long double x) 162 long double log1pl(long double x) {
174 { 163 return log1p(x);
175 » return log1p(x);
176 } 164 }
177 #endif 165 #endif
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698