1 /* log10l.c 2 * 3 * Common logarithm, 128-bit long double precision 4 * 5 * 6 * 7 * SYNOPSIS: 8 * 9 * long double x, y, log10l(); 10 * 11 * y = log10l( x ); 12 * 13 * 14 * 15 * DESCRIPTION: 16 * 17 * Returns the base 10 logarithm of x. 18 * 19 * The argument is separated into its exponent and fractional 20 * parts. If the exponent is between -1 and +1, the logarithm 21 * of the fraction is approximated by 22 * 23 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x). 24 * 25 * Otherwise, setting z = 2(x-1)/x+1), 26 * 27 * log(x) = z + z^3 P(z)/Q(z). 28 * 29 * 30 * 31 * ACCURACY: 32 * 33 * Relative error: 34 * arithmetic domain # trials peak rms 35 * IEEE 0.5, 2.0 30000 2.3e-34 4.9e-35 36 * IEEE exp(+-10000) 30000 1.0e-34 4.1e-35 37 * 38 * In the tests over the interval exp(+-10000), the logarithms 39 * of the random arguments were uniformly distributed over 40 * [-10000, +10000]. 41 * 42 */ 43 44 /* 45 Cephes Math Library Release 2.2: January, 1991 46 Copyright 1984, 1991 by Stephen L. Moshier 47 Adapted for glibc November, 2001 48 49 This library is free software; you can redistribute it and/or 50 modify it under the terms of the GNU Lesser General Public 51 License as published by the Free Software Foundation; either 52 version 2.1 of the License, or (at your option) any later version. 53 54 This library is distributed in the hope that it will be useful, 55 but WITHOUT ANY WARRANTY; without even the implied warranty of 56 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 57 Lesser General Public License for more details. 58 59 You should have received a copy of the GNU Lesser General Public 60 License along with this library; if not, see <http://www.gnu.org/licenses/>. 61 */ 62 63 #include "quadmath-imp.h" 64 65 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x) 66 * 1/sqrt(2) <= x < sqrt(2) 67 * Theoretical peak relative error = 5.3e-37, 68 * relative peak error spread = 2.3e-14 69 */ 70 static const __float128 P[13] = 71 { 72 1.313572404063446165910279910527789794488E4Q, 73 7.771154681358524243729929227226708890930E4Q, 74 2.014652742082537582487669938141683759923E5Q, 75 3.007007295140399532324943111654767187848E5Q, 76 2.854829159639697837788887080758954924001E5Q, 77 1.797628303815655343403735250238293741397E5Q, 78 7.594356839258970405033155585486712125861E4Q, 79 2.128857716871515081352991964243375186031E4Q, 80 3.824952356185897735160588078446136783779E3Q, 81 4.114517881637811823002128927449878962058E2Q, 82 2.321125933898420063925789532045674660756E1Q, 83 4.998469661968096229986658302195402690910E-1Q, 84 1.538612243596254322971797716843006400388E-6Q 85 }; 86 static const __float128 Q[12] = 87 { 88 3.940717212190338497730839731583397586124E4Q, 89 2.626900195321832660448791748036714883242E5Q, 90 7.777690340007566932935753241556479363645E5Q, 91 1.347518538384329112529391120390701166528E6Q, 92 1.514882452993549494932585972882995548426E6Q, 93 1.158019977462989115839826904108208787040E6Q, 94 6.132189329546557743179177159925690841200E5Q, 95 2.248234257620569139969141618556349415120E5Q, 96 5.605842085972455027590989944010492125825E4Q, 97 9.147150349299596453976674231612674085381E3Q, 98 9.104928120962988414618126155557301584078E2Q, 99 4.839208193348159620282142911143429644326E1Q 100 /* 1.000000000000000000000000000000000000000E0L, */ 101 }; 102 103 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), 104 * where z = 2(x-1)/(x+1) 105 * 1/sqrt(2) <= x < sqrt(2) 106 * Theoretical peak relative error = 1.1e-35, 107 * relative peak error spread 1.1e-9 108 */ 109 static const __float128 R[6] = 110 { 111 1.418134209872192732479751274970992665513E5Q, 112 -8.977257995689735303686582344659576526998E4Q, 113 2.048819892795278657810231591630928516206E4Q, 114 -2.024301798136027039250415126250455056397E3Q, 115 8.057002716646055371965756206836056074715E1Q, 116 -8.828896441624934385266096344596648080902E-1Q 117 }; 118 static const __float128 S[6] = 119 { 120 1.701761051846631278975701529965589676574E6Q, 121 -1.332535117259762928288745111081235577029E6Q, 122 4.001557694070773974936904547424676279307E5Q, 123 -5.748542087379434595104154610899551484314E4Q, 124 3.998526750980007367835804959888064681098E3Q, 125 -1.186359407982897997337150403816839480438E2Q 126 /* 1.000000000000000000000000000000000000000E0L, */ 127 }; 128 129 static const __float128 130 /* log10(2) */ 131 L102A = 0.3125Q, 132 L102B = -1.14700043360188047862611052755069732318101185E-2Q, 133 /* log10(e) */ 134 L10EA = 0.5Q, 135 L10EB = -6.570551809674817234887108108339491770560299E-2Q, 136 /* sqrt(2)/2 */ 137 SQRTH = 7.071067811865475244008443621048490392848359E-1Q; 138 139 140 141 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ 142 143 static __float128 144 neval (__float128 x, const __float128 *p, int n) 145 { 146 __float128 y; 147 148 p += n; 149 y = *p--; 150 do 151 { 152 y = y * x + *p--; 153 } 154 while (--n > 0); 155 return y; 156 } 157 158 159 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */ 160 161 static __float128 162 deval (__float128 x, const __float128 *p, int n) 163 { 164 __float128 y; 165 166 p += n; 167 y = x + *p--; 168 do 169 { 170 y = y * x + *p--; 171 } 172 while (--n > 0); 173 return y; 174 } 175 176 177 178 __float128 179 log10q (__float128 x) 180 { 181 __float128 z; 182 __float128 y; 183 int e; 184 int64_t hx, lx; 185 186 /* Test for domain */ 187 GET_FLT128_WORDS64 (hx, lx, x); 188 if (((hx & 0x7fffffffffffffffLL) | lx) == 0) 189 return (-1 / fabsq (x)); /* log10l(+-0)=-inf */ 190 if (hx < 0) 191 return (x - x) / (x - x); 192 if (hx >= 0x7fff000000000000LL) 193 return (x + x); 194 195 if (x == 1) 196 return 0; 197 198 /* separate mantissa from exponent */ 199 200 /* Note, frexp is used so that denormal numbers 201 * will be handled properly. 202 */ 203 x = frexpq (x, &e); 204 205 206 /* logarithm using log(x) = z + z**3 P(z)/Q(z), 207 * where z = 2(x-1)/x+1) 208 */ 209 if ((e > 2) || (e < -2)) 210 { 211 if (x < SQRTH) 212 { /* 2( 2x-1 )/( 2x+1 ) */ 213 e -= 1; 214 z = x - 0.5Q; 215 y = 0.5Q * z + 0.5Q; 216 } 217 else 218 { /* 2 (x-1)/(x+1) */ 219 z = x - 0.5Q; 220 z -= 0.5Q; 221 y = 0.5Q * x + 0.5Q; 222 } 223 x = z / y; 224 z = x * x; 225 y = x * (z * neval (z, R, 5) / deval (z, S, 5)); 226 goto done; 227 } 228 229 230 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ 231 232 if (x < SQRTH) 233 { 234 e -= 1; 235 x = 2.0 * x - 1; /* 2x - 1 */ 236 } 237 else 238 { 239 x = x - 1; 240 } 241 z = x * x; 242 y = x * (z * neval (x, P, 12) / deval (x, Q, 11)); 243 y = y - 0.5 * z; 244 245 done: 246 247 /* Multiply log of fraction by log10(e) 248 * and base 2 exponent by log10(2). 249 */ 250 z = y * L10EB; 251 z += x * L10EB; 252 z += e * L102B; 253 z += y * L10EA; 254 z += x * L10EA; 255 z += e * L102A; 256 return (z); 257 } 258