1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_logl.c */ 2 /* 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 4 * 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 7 * copyright notice and this permission notice appear in all copies. 8 * 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 11 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 12 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 13 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 14 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 15 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 16 */ 17 /* 18 * Natural logarithm, long double precision 19 * 20 * 21 * SYNOPSIS: 22 * 23 * long double x, y, logl(); 24 * 25 * y = logl( x ); 26 * 27 * 28 * DESCRIPTION: 29 * 30 * Returns the base e (2.718...) logarithm of x. 31 * 32 * The argument is separated into its exponent and fractional 33 * parts. If the exponent is between -1 and +1, the logarithm 34 * of the fraction is approximated by 35 * 36 * log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x). 37 * 38 * Otherwise, setting z = 2(x-1)/(x+1), 39 * 40 * log(x) = log(1+z/2) - log(1-z/2) = z + z**3 P(z)/Q(z). 41 * 42 * 43 * ACCURACY: 44 * 45 * Relative error: 46 * arithmetic domain # trials peak rms 47 * IEEE 0.5, 2.0 150000 8.71e-20 2.75e-20 48 * IEEE exp(+-10000) 100000 5.39e-20 2.34e-20 49 * 50 * In the tests over the interval exp(+-10000), the logarithms 51 * of the random arguments were uniformly distributed over 52 * [-10000, +10000]. 53 */ 54 55 #include "libm.h" 56 57 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 58 long double logl(long double x) 59 { 60 return log(x); 61 } 62 #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) 64 * 1/sqrt(2) <= x < sqrt(2) 65 * Theoretical peak relative error = 2.32e-20 66 */ 67 static const long double P[] = { 68 4.5270000862445199635215E-5L, 69 4.9854102823193375972212E-1L, 70 6.5787325942061044846969E0L, 71 2.9911919328553073277375E1L, 72 6.0949667980987787057556E1L, 73 5.7112963590585538103336E1L, 74 2.0039553499201281259648E1L, 75 }; 76 static const long double Q[] = { 77 /* 1.0000000000000000000000E0,*/ 78 1.5062909083469192043167E1L, 79 8.3047565967967209469434E1L, 80 2.2176239823732856465394E2L, 81 3.0909872225312059774938E2L, 82 2.1642788614495947685003E2L, 83 6.0118660497603843919306E1L, 84 }; 85 86 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), 87 * where z = 2(x-1)/(x+1) 88 * 1/sqrt(2) <= x < sqrt(2) 89 * Theoretical peak relative error = 6.16e-22 90 */ 91 static const long double R[4] = { 92 1.9757429581415468984296E-3L, 93 -7.1990767473014147232598E-1L, 94 1.0777257190312272158094E1L, 95 -3.5717684488096787370998E1L, 96 }; 97 static const long double S[4] = { 98 /* 1.00000000000000000000E0L,*/ 99 -2.6201045551331104417768E1L, 100 1.9361891836232102174846E2L, 101 -4.2861221385716144629696E2L, 102 }; 103 static const long double C1 = 6.9314575195312500000000E-1L; 104 static const long double C2 = 1.4286068203094172321215E-6L; 105 106 #define SQRTH 0.70710678118654752440L 107 108 long double logl(long double x) 109 { 110 long double y, z; 111 int e; 112 113 if (isnan(x)) 114 return x; 115 if (x == INFINITY) 116 return x; 117 if (x <= 0.0) { 118 if (x == 0.0) 119 return -1/(x*x); /* -inf with divbyzero */ 120 return 0/0.0f; /* nan with invalid */ 121 } 122 123 /* separate mantissa from exponent */ 124 /* Note, frexp is used so that denormal numbers 125 * will be handled properly. 126 */ 127 x = frexpl(x, &e); 128 129 /* logarithm using log(x) = z + z**3 P(z)/Q(z), 130 * where z = 2(x-1)/(x+1) 131 */ 132 if (e > 2 || e < -2) { 133 if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ 134 e -= 1; 135 z = x - 0.5; 136 y = 0.5 * z + 0.5; 137 } else { /* 2 (x-1)/(x+1) */ 138 z = x - 0.5; 139 z -= 0.5; 140 y = 0.5 * x + 0.5; 141 } 142 x = z / y; 143 z = x*x; 144 z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); 145 z = z + e * C2; 146 z = z + x; 147 z = z + e * C1; 148 return z; 149 } 150 151 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ 152 if (x < SQRTH) { 153 e -= 1; 154 x = 2.0*x - 1.0; 155 } else { 156 x = x - 1.0; 157 } 158 z = x*x; 159 y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6)); 160 y = y + e * C2; 161 z = y - 0.5*z; 162 /* Note, the sum of above terms does not exceed x/4, 163 * so it contributes at most about 1/4 lsb to the error. 164 */ 165 z = z + x; 166 z = z + e * C1; /* This sum has an error of 1/2 lsb. */ 167 return z; 168 } 169 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 170 // TODO: broken implementation to make things compile 171 long double logl(long double x) 172 { 173 return log(x); 174 } 175 #endif 176