OLD | NEW |
1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_logl.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_logl.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 37 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
48 * IEEE exp(+-10000) 100000 5.39e-20 2.34e-20 | 48 * IEEE exp(+-10000) 100000 5.39e-20 2.34e-20 |
49 * | 49 * |
50 * In the tests over the interval exp(+-10000), the logarithms | 50 * In the tests over the interval exp(+-10000), the logarithms |
51 * of the random arguments were uniformly distributed over | 51 * of the random arguments were uniformly distributed over |
52 * [-10000, +10000]. | 52 * [-10000, +10000]. |
53 */ | 53 */ |
54 | 54 |
55 #include "libm.h" | 55 #include "libm.h" |
56 | 56 |
57 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 57 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
58 long double logl(long double x) | 58 long double logl(long double x) { |
59 { | 59 return log(x); |
60 » return log(x); | |
61 } | 60 } |
62 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 | 61 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 |
63 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x) | 62 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x) |
64 * 1/sqrt(2) <= x < sqrt(2) | 63 * 1/sqrt(2) <= x < sqrt(2) |
65 * Theoretical peak relative error = 2.32e-20 | 64 * Theoretical peak relative error = 2.32e-20 |
66 */ | 65 */ |
67 static const long double P[] = { | 66 static const long double P[] = { |
68 4.5270000862445199635215E-5L, | 67 4.5270000862445199635215E-5L, 4.9854102823193375972212E-1L, |
69 4.9854102823193375972212E-1L, | 68 6.5787325942061044846969E0L, 2.9911919328553073277375E1L, |
70 6.5787325942061044846969E0L, | 69 6.0949667980987787057556E1L, 5.7112963590585538103336E1L, |
71 2.9911919328553073277375E1L, | 70 2.0039553499201281259648E1L, |
72 6.0949667980987787057556E1L, | |
73 5.7112963590585538103336E1L, | |
74 2.0039553499201281259648E1L, | |
75 }; | 71 }; |
76 static const long double Q[] = { | 72 static const long double Q[] = { |
77 /* 1.0000000000000000000000E0,*/ | 73 /* 1.0000000000000000000000E0,*/ |
78 1.5062909083469192043167E1L, | 74 1.5062909083469192043167E1L, 8.3047565967967209469434E1L, |
79 8.3047565967967209469434E1L, | 75 2.2176239823732856465394E2L, 3.0909872225312059774938E2L, |
80 2.2176239823732856465394E2L, | 76 2.1642788614495947685003E2L, 6.0118660497603843919306E1L, |
81 3.0909872225312059774938E2L, | |
82 2.1642788614495947685003E2L, | |
83 6.0118660497603843919306E1L, | |
84 }; | 77 }; |
85 | 78 |
86 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), | 79 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), |
87 * where z = 2(x-1)/(x+1) | 80 * where z = 2(x-1)/(x+1) |
88 * 1/sqrt(2) <= x < sqrt(2) | 81 * 1/sqrt(2) <= x < sqrt(2) |
89 * Theoretical peak relative error = 6.16e-22 | 82 * Theoretical peak relative error = 6.16e-22 |
90 */ | 83 */ |
91 static const long double R[4] = { | 84 static const long double R[4] = { |
92 1.9757429581415468984296E-3L, | 85 1.9757429581415468984296E-3L, -7.1990767473014147232598E-1L, |
93 -7.1990767473014147232598E-1L, | 86 1.0777257190312272158094E1L, -3.5717684488096787370998E1L, |
94 1.0777257190312272158094E1L, | |
95 -3.5717684488096787370998E1L, | |
96 }; | 87 }; |
97 static const long double S[4] = { | 88 static const long double S[4] = { |
98 /* 1.00000000000000000000E0L,*/ | 89 /* 1.00000000000000000000E0L,*/ |
99 -2.6201045551331104417768E1L, | 90 -2.6201045551331104417768E1L, 1.9361891836232102174846E2L, |
100 1.9361891836232102174846E2L, | 91 -4.2861221385716144629696E2L, |
101 -4.2861221385716144629696E2L, | |
102 }; | 92 }; |
103 static const long double C1 = 6.9314575195312500000000E-1L; | 93 static const long double C1 = 6.9314575195312500000000E-1L; |
104 static const long double C2 = 1.4286068203094172321215E-6L; | 94 static const long double C2 = 1.4286068203094172321215E-6L; |
105 | 95 |
106 #define SQRTH 0.70710678118654752440L | 96 #define SQRTH 0.70710678118654752440L |
107 | 97 |
108 long double logl(long double x) | 98 long double logl(long double x) { |
109 { | 99 long double y, z; |
110 » long double y, z; | 100 int e; |
111 » int e; | |
112 | 101 |
113 » if (isnan(x)) | 102 if (isnan(x)) |
114 » » return x; | 103 return x; |
115 » if (x == INFINITY) | 104 if (x == INFINITY) |
116 » » return x; | 105 return x; |
117 » if (x <= 0.0) { | 106 if (x <= 0.0) { |
118 » » if (x == 0.0) | 107 if (x == 0.0) |
119 » » » return -1/(x*x); /* -inf with divbyzero */ | 108 return -1 / (x * x); /* -inf with divbyzero */ |
120 » » return 0/0.0f; /* nan with invalid */ | 109 return 0 / 0.0f; /* nan with invalid */ |
121 » } | 110 } |
122 | 111 |
123 » /* separate mantissa from exponent */ | 112 /* separate mantissa from exponent */ |
124 » /* Note, frexp is used so that denormal numbers | 113 /* Note, frexp is used so that denormal numbers |
125 » * will be handled properly. | 114 * will be handled properly. |
126 » */ | 115 */ |
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 » */ | 120 */ |
132 » if (e > 2 || e < -2) { | 121 if (e > 2 || e < -2) { |
133 » » if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ | 122 if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ |
134 » » » e -= 1; | 123 e -= 1; |
135 » » » z = x - 0.5; | 124 z = x - 0.5; |
136 » » » y = 0.5 * z + 0.5; | 125 y = 0.5 * z + 0.5; |
137 » » } else { /* 2 (x-1)/(x+1) */ | 126 } else { /* 2 (x-1)/(x+1) */ |
138 » » » z = x - 0.5; | 127 z = x - 0.5; |
139 » » » z -= 0.5; | 128 z -= 0.5; |
140 » » » y = 0.5 * x + 0.5; | 129 y = 0.5 * x + 0.5; |
141 » » } | 130 } |
142 » » x = z / y; | 131 x = z / y; |
143 » » z = x*x; | 132 z = x * x; |
144 » » z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); | 133 z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); |
145 » » z = z + e * C2; | 134 z = z + e * C2; |
146 » » z = z + x; | 135 z = z + x; |
147 » » z = z + e * C1; | 136 z = z + e * C1; |
148 » » return z; | 137 return z; |
149 » } | 138 } |
150 | 139 |
151 » /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ | 140 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ |
152 » if (x < SQRTH) { | 141 if (x < SQRTH) { |
153 » » e -= 1; | 142 e -= 1; |
154 » » x = 2.0*x - 1.0; | 143 x = 2.0 * x - 1.0; |
155 » } else { | 144 } else { |
156 » » x = x - 1.0; | 145 x = x - 1.0; |
157 » } | 146 } |
158 » z = x*x; | 147 z = x * x; |
159 » y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6)); | 148 y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6)); |
160 » y = y + e * C2; | 149 y = y + e * C2; |
161 » z = y - 0.5*z; | 150 z = y - 0.5 * z; |
162 » /* Note, the sum of above terms does not exceed x/4, | 151 /* Note, the sum of above terms does not exceed x/4, |
163 » * so it contributes at most about 1/4 lsb to the error. | 152 * so it contributes at most about 1/4 lsb to the error. |
164 » */ | 153 */ |
165 » z = z + x; | 154 z = z + x; |
166 » z = z + e * C1; /* This sum has an error of 1/2 lsb. */ | 155 z = z + e * C1; /* This sum has an error of 1/2 lsb. */ |
167 » return z; | 156 return z; |
168 } | 157 } |
169 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 | 158 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 |
170 // TODO: broken implementation to make things compile | 159 // TODO: broken implementation to make things compile |
171 long double logl(long double x) | 160 long double logl(long double x) { |
172 { | 161 return log(x); |
173 » return log(x); | |
174 } | 162 } |
175 #endif | 163 #endif |
OLD | NEW |