1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_log1pl.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 * Relative error logarithm 19 * Natural logarithm of 1+x, long double precision 20 * 21 * 22 * SYNOPSIS: 23 * 24 * long double x, y, log1pl(); 25 * 26 * y = log1pl( x ); 27 * 28 * 29 * DESCRIPTION: 30 * 31 * Returns the base e (2.718...) logarithm of 1+x. 32 * 33 * The argument 1+x is separated into its exponent and fractional 34 * parts. If the exponent is between -1 and +1, the logarithm 35 * of the fraction is approximated by 36 * 37 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x). 38 * 39 * Otherwise, setting z = 2(x-1)/x+1), 40 * 41 * log(x) = z + z^3 P(z)/Q(z). 42 * 43 * 44 * ACCURACY: 45 * 46 * Relative error: 47 * arithmetic domain # trials peak rms 48 * IEEE -1.0, 9.0 100000 8.2e-20 2.5e-20 49 */ 50 51 #include "libm.h" 52 53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 54 long double log1pl(long double x) 55 { 56 return log1p(x); 57 } 58 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 59 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x) 60 * 1/sqrt(2) <= x < sqrt(2) 61 * Theoretical peak relative error = 2.32e-20 62 */ 63 static const long double P[] = { 64 4.5270000862445199635215E-5L, 65 4.9854102823193375972212E-1L, 66 6.5787325942061044846969E0L, 67 2.9911919328553073277375E1L, 68 6.0949667980987787057556E1L, 69 5.7112963590585538103336E1L, 70 2.0039553499201281259648E1L, 71 }; 72 static const long double Q[] = { 73 /* 1.0000000000000000000000E0,*/ 74 1.5062909083469192043167E1L, 75 8.3047565967967209469434E1L, 76 2.2176239823732856465394E2L, 77 3.0909872225312059774938E2L, 78 2.1642788614495947685003E2L, 79 6.0118660497603843919306E1L, 80 }; 81 82 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), 83 * where z = 2(x-1)/(x+1) 84 * 1/sqrt(2) <= x < sqrt(2) 85 * Theoretical peak relative error = 6.16e-22 86 */ 87 static const long double R[4] = { 88 1.9757429581415468984296E-3L, 89 -7.1990767473014147232598E-1L, 90 1.0777257190312272158094E1L, 91 -3.5717684488096787370998E1L, 92 }; 93 static const long double S[4] = { 94 /* 1.00000000000000000000E0L,*/ 95 -2.6201045551331104417768E1L, 96 1.9361891836232102174846E2L, 97 -4.2861221385716144629696E2L, 98 }; 99 static const long double C1 = 6.9314575195312500000000E-1L; 100 static const long double C2 = 1.4286068203094172321215E-6L; 101 102 #define SQRTH 0.70710678118654752440L 103 104 long double log1pl(long double xm1) 105 { 106 long double x, y, z; 107 int e; 108 109 if (isnan(xm1)) 110 return xm1; 111 if (xm1 == INFINITY) 112 return xm1; 113 if (xm1 == 0.0) 114 return xm1; 115 116 x = xm1 + 1.0; 117 118 /* Test for domain errors. */ 119 if (x <= 0.0) { 120 if (x == 0.0) 121 return -1/(x*x); /* -inf with divbyzero */ 122 return 0/0.0f; /* nan with invalid */ 123 } 124 125 /* Separate mantissa from exponent. 126 Use frexp so that denormal numbers will be handled properly. */ 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 if (e > 2 || e < -2) { 132 if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ 133 e -= 1; 134 z = x - 0.5; 135 y = 0.5 * z + 0.5; 136 } else { /* 2 (x-1)/(x+1) */ 137 z = x - 0.5; 138 z -= 0.5; 139 y = 0.5 * x + 0.5; 140 } 141 x = z / y; 142 z = x*x; 143 z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); 144 z = z + e * C2; 145 z = z + x; 146 z = z + e * C1; 147 return z; 148 } 149 150 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ 151 if (x < SQRTH) { 152 e -= 1; 153 if (e != 0) 154 x = 2.0 * x - 1.0; 155 else 156 x = xm1; 157 } else { 158 if (e != 0) 159 x = x - 1.0; 160 else 161 x = xm1; 162 } 163 z = x*x; 164 y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6)); 165 y = y + e * C2; 166 z = y - 0.5 * z; 167 z = z + x; 168 z = z + e * C1; 169 return z; 170 } 171 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 172 // TODO: broken implementation to make things compile 173 long double log1pl(long double x) 174 { 175 return log1p(x); 176 } 177 #endif 178