1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_log2l.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 * Base 2 logarithm, long double precision 19 * 20 * 21 * SYNOPSIS: 22 * 23 * long double x, y, log2l(); 24 * 25 * y = log2l( x ); 26 * 27 * 28 * DESCRIPTION: 29 * 30 * Returns the base 2 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 (natural) 34 * logarithm 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) = 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 30000 9.8e-20 2.7e-20 48 * IEEE exp(+-10000) 70000 5.4e-20 2.3e-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 log2l(long double x) 59 { 60 return log2(x); 61 } 62 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 63 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x) 64 * 1/sqrt(2) <= x < sqrt(2) 65 * Theoretical peak relative error = 6.2e-22 66 */ 67 static const long double P[] = { 68 4.9962495940332550844739E-1L, 69 1.0767376367209449010438E1L, 70 7.7671073698359539859595E1L, 71 2.5620629828144409632571E2L, 72 4.2401812743503691187826E2L, 73 3.4258224542413922935104E2L, 74 1.0747524399916215149070E2L, 75 }; 76 static const long double Q[] = { 77 /* 1.0000000000000000000000E0,*/ 78 2.3479774160285863271658E1L, 79 1.9444210022760132894510E2L, 80 7.7952888181207260646090E2L, 81 1.6911722418503949084863E3L, 82 2.0307734695595183428202E3L, 83 1.2695660352705325274404E3L, 84 3.2242573199748645407652E2L, 85 }; 86 87 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), 88 * where z = 2(x-1)/(x+1) 89 * 1/sqrt(2) <= x < sqrt(2) 90 * Theoretical peak relative error = 6.16e-22 91 */ 92 static const long double R[4] = { 93 1.9757429581415468984296E-3L, 94 -7.1990767473014147232598E-1L, 95 1.0777257190312272158094E1L, 96 -3.5717684488096787370998E1L, 97 }; 98 static const long double S[4] = { 99 /* 1.00000000000000000000E0L,*/ 100 -2.6201045551331104417768E1L, 101 1.9361891836232102174846E2L, 102 -4.2861221385716144629696E2L, 103 }; 104 /* log2(e) - 1 */ 105 #define LOG2EA 4.4269504088896340735992e-1L 106 107 #define SQRTH 0.70710678118654752440L 108 109 long double log2l(long double x) 110 { 111 long double y, z; 112 int e; 113 114 if (isnan(x)) 115 return x; 116 if (x == INFINITY) 117 return x; 118 if (x <= 0.0) { 119 if (x == 0.0) 120 return -1/(x*x); /* -inf with divbyzero */ 121 return 0/0.0f; /* nan with invalid */ 122 } 123 124 /* separate mantissa from exponent */ 125 /* Note, frexp is used so that denormal numbers 126 * will be handled properly. 127 */ 128 x = frexpl(x, &e); 129 130 /* logarithm using log(x) = z + z**3 P(z)/Q(z), 131 * where z = 2(x-1)/x+1) 132 */ 133 if (e > 2 || e < -2) { 134 if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ 135 e -= 1; 136 z = x - 0.5; 137 y = 0.5 * z + 0.5; 138 } else { /* 2 (x-1)/(x+1) */ 139 z = x - 0.5; 140 z -= 0.5; 141 y = 0.5 * x + 0.5; 142 } 143 x = z / y; 144 z = x*x; 145 y = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); 146 goto done; 147 } 148 149 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ 150 if (x < SQRTH) { 151 e -= 1; 152 x = 2.0*x - 1.0; 153 } else { 154 x = x - 1.0; 155 } 156 z = x*x; 157 y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7)); 158 y = y - 0.5*z; 159 160 done: 161 /* Multiply log of fraction by log2(e) 162 * and base 2 exponent by 1 163 * 164 * ***CAUTION*** 165 * 166 * This sequence of operations is critical and it may 167 * be horribly defeated by some compiler optimizers. 168 */ 169 z = y * LOG2EA; 170 z += x * LOG2EA; 171 z += y; 172 z += x; 173 z += e; 174 return z; 175 } 176 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 177 // TODO: broken implementation to make things compile 178 long double log2l(long double x) 179 { 180 return log2(x); 181 } 182 #endif 183