184447Sdfr/* log1pq.c 284447Sdfr * 384447Sdfr * Relative error logarithm 484447Sdfr * Natural logarithm of 1+x, 128-bit long double precision 584447Sdfr * 684447Sdfr * 784447Sdfr * 884447Sdfr * SYNOPSIS: 984447Sdfr * 1084447Sdfr * long double x, y, log1pq(); 1184447Sdfr * 1284447Sdfr * y = log1pq( x ); 1384447Sdfr * 1484447Sdfr * 1584447Sdfr * 1684447Sdfr * DESCRIPTION: 1784447Sdfr * 1884447Sdfr * Returns the base e (2.718...) logarithm of 1+x. 1984447Sdfr * 2084447Sdfr * The argument 1+x is separated into its exponent and fractional 2184447Sdfr * parts. If the exponent is between -1 and +1, the logarithm 2284447Sdfr * of the fraction is approximated by 2384447Sdfr * 2484447Sdfr * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x). 2584447Sdfr * 2684447Sdfr * Otherwise, setting z = 2(w-1)/(w+1), 2784447Sdfr * 2884447Sdfr * log(w) = z + z^3 P(z)/Q(z). 2984447Sdfr * 3084447Sdfr * 3184447Sdfr * 32193530Sjkim * ACCURACY: 33193530Sjkim * 3484447Sdfr * Relative error: 3584447Sdfr * arithmetic domain # trials peak rms 3684447Sdfr * IEEE -1, 8 100000 1.9e-34 4.3e-35 3784447Sdfr */ 3884447Sdfr 3984447Sdfr/* Copyright 2001 by Stephen L. Moshier 4084447Sdfr 4184447Sdfr This library is free software; you can redistribute it and/or 4284447Sdfr modify it under the terms of the GNU Lesser General Public 4384447Sdfr License as published by the Free Software Foundation; either 4484447Sdfr version 2.1 of the License, or (at your option) any later version. 4584447Sdfr 46 This library is distributed in the hope that it will be useful, 47 but WITHOUT ANY WARRANTY; without even the implied warranty of 48 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 49 Lesser General Public License for more details. 50 51 You should have received a copy of the GNU Lesser General Public 52 License along with this library; if not, see 53 <http://www.gnu.org/licenses/>. */ 54 55#include "quadmath-imp.h" 56 57/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x) 58 * 1/sqrt(2) <= 1+x < sqrt(2) 59 * Theoretical peak relative error = 5.3e-37, 60 * relative peak error spread = 2.3e-14 61 */ 62static const __float128 63 P12 = 1.538612243596254322971797716843006400388E-6Q, 64 P11 = 4.998469661968096229986658302195402690910E-1Q, 65 P10 = 2.321125933898420063925789532045674660756E1Q, 66 P9 = 4.114517881637811823002128927449878962058E2Q, 67 P8 = 3.824952356185897735160588078446136783779E3Q, 68 P7 = 2.128857716871515081352991964243375186031E4Q, 69 P6 = 7.594356839258970405033155585486712125861E4Q, 70 P5 = 1.797628303815655343403735250238293741397E5Q, 71 P4 = 2.854829159639697837788887080758954924001E5Q, 72 P3 = 3.007007295140399532324943111654767187848E5Q, 73 P2 = 2.014652742082537582487669938141683759923E5Q, 74 P1 = 7.771154681358524243729929227226708890930E4Q, 75 P0 = 1.313572404063446165910279910527789794488E4Q, 76 /* Q12 = 1.000000000000000000000000000000000000000E0L, */ 77 Q11 = 4.839208193348159620282142911143429644326E1Q, 78 Q10 = 9.104928120962988414618126155557301584078E2Q, 79 Q9 = 9.147150349299596453976674231612674085381E3Q, 80 Q8 = 5.605842085972455027590989944010492125825E4Q, 81 Q7 = 2.248234257620569139969141618556349415120E5Q, 82 Q6 = 6.132189329546557743179177159925690841200E5Q, 83 Q5 = 1.158019977462989115839826904108208787040E6Q, 84 Q4 = 1.514882452993549494932585972882995548426E6Q, 85 Q3 = 1.347518538384329112529391120390701166528E6Q, 86 Q2 = 7.777690340007566932935753241556479363645E5Q, 87 Q1 = 2.626900195321832660448791748036714883242E5Q, 88 Q0 = 3.940717212190338497730839731583397586124E4Q; 89 90/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), 91 * where z = 2(x-1)/(x+1) 92 * 1/sqrt(2) <= x < sqrt(2) 93 * Theoretical peak relative error = 1.1e-35, 94 * relative peak error spread 1.1e-9 95 */ 96static const __float128 97 R5 = -8.828896441624934385266096344596648080902E-1Q, 98 R4 = 8.057002716646055371965756206836056074715E1Q, 99 R3 = -2.024301798136027039250415126250455056397E3Q, 100 R2 = 2.048819892795278657810231591630928516206E4Q, 101 R1 = -8.977257995689735303686582344659576526998E4Q, 102 R0 = 1.418134209872192732479751274970992665513E5Q, 103 /* S6 = 1.000000000000000000000000000000000000000E0L, */ 104 S5 = -1.186359407982897997337150403816839480438E2Q, 105 S4 = 3.998526750980007367835804959888064681098E3Q, 106 S3 = -5.748542087379434595104154610899551484314E4Q, 107 S2 = 4.001557694070773974936904547424676279307E5Q, 108 S1 = -1.332535117259762928288745111081235577029E6Q, 109 S0 = 1.701761051846631278975701529965589676574E6Q; 110 111/* C1 + C2 = ln 2 */ 112static const __float128 C1 = 6.93145751953125E-1Q; 113static const __float128 C2 = 1.428606820309417232121458176568075500134E-6Q; 114 115static const __float128 sqrth = 0.7071067811865475244008443621048490392848Q; 116/* ln (2^16384 * (1 - 2^-113)) */ 117static const __float128 zero = 0; 118 119__float128 120log1pq (__float128 xm1) 121{ 122 __float128 x, y, z, r, s; 123 ieee854_float128 u; 124 int32_t hx; 125 int e; 126 127 /* Test for NaN or infinity input. */ 128 u.value = xm1; 129 hx = u.words32.w0; 130 if ((hx & 0x7fffffff) >= 0x7fff0000) 131 return xm1 + fabsq (xm1); 132 133 /* log1p(+- 0) = +- 0. */ 134 if (((hx & 0x7fffffff) == 0) 135 && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) 136 return xm1; 137 138 if ((hx & 0x7fffffff) < 0x3f8e0000) 139 { 140 math_check_force_underflow (xm1); 141 if ((int) xm1 == 0) 142 return xm1; 143 } 144 145 if (xm1 >= 0x1p113Q) 146 x = xm1; 147 else 148 x = xm1 + 1; 149 150 /* log1p(-1) = -inf */ 151 if (x <= 0) 152 { 153 if (x == 0) 154 return (-1 / zero); /* log1p(-1) = -inf */ 155 else 156 return (zero / (x - x)); 157 } 158 159 /* Separate mantissa from exponent. */ 160 161 /* Use frexp used so that denormal numbers will be handled properly. */ 162 x = frexpq (x, &e); 163 164 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2), 165 where z = 2(x-1)/x+1). */ 166 if ((e > 2) || (e < -2)) 167 { 168 if (x < sqrth) 169 { /* 2( 2x-1 )/( 2x+1 ) */ 170 e -= 1; 171 z = x - 0.5Q; 172 y = 0.5Q * z + 0.5Q; 173 } 174 else 175 { /* 2 (x-1)/(x+1) */ 176 z = x - 0.5Q; 177 z -= 0.5Q; 178 y = 0.5Q * x + 0.5Q; 179 } 180 x = z / y; 181 z = x * x; 182 r = ((((R5 * z 183 + R4) * z 184 + R3) * z 185 + R2) * z 186 + R1) * z 187 + R0; 188 s = (((((z 189 + S5) * z 190 + S4) * z 191 + S3) * z 192 + S2) * z 193 + S1) * z 194 + S0; 195 z = x * (z * r / s); 196 z = z + e * C2; 197 z = z + x; 198 z = z + e * C1; 199 return (z); 200 } 201 202 203 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */ 204 205 if (x < sqrth) 206 { 207 e -= 1; 208 if (e != 0) 209 x = 2 * x - 1; /* 2x - 1 */ 210 else 211 x = xm1; 212 } 213 else 214 { 215 if (e != 0) 216 x = x - 1; 217 else 218 x = xm1; 219 } 220 z = x * x; 221 r = (((((((((((P12 * x 222 + P11) * x 223 + P10) * x 224 + P9) * x 225 + P8) * x 226 + P7) * x 227 + P6) * x 228 + P5) * x 229 + P4) * x 230 + P3) * x 231 + P2) * x 232 + P1) * x 233 + P0; 234 s = (((((((((((x 235 + Q11) * x 236 + Q10) * x 237 + Q9) * x 238 + Q8) * x 239 + Q7) * x 240 + Q6) * x 241 + Q5) * x 242 + Q4) * x 243 + Q3) * x 244 + Q2) * x 245 + Q1) * x 246 + Q0; 247 y = x * (z * r / s); 248 y = y + e * C2; 249 z = y - 0.5Q * z; 250 z = z + x; 251 z = z + e * C1; 252 return (z); 253} 254