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