OLD | NEW |
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 Loading... |
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 |
OLD | NEW |