1 /* 2 * Copyright (c) 1992 Regents of the University of California. 3 * All rights reserved. 4 * 5 * %sccs.include.redist.c% 6 */ 7 8 #ifndef lint 9 static char sccsid[] = "@(#)log.c 5.9 (Berkeley) 12/16/92"; 10 #endif /* not lint */ 11 12 #include <math.h> 13 #include <errno.h> 14 15 #include "log_table.h" 16 #include "mathimpl.h" 17 18 /* Table-driven natural logarithm. 19 * 20 * This code was derived, with minor modifications, from: 21 * Peter Tang, "Table-Driven Implementation of the 22 * Logarithm in IEEE Floating-Point arithmetic." ACM Trans. 23 * Math Software, vol 16. no 4, pp 378-400, Dec 1990). 24 * 25 * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256, 26 * where F = j/128 for j an integer in [0, 128]. 27 * 28 * log(2^m) = log2_hi*m + log2_tail*m 29 * since m is an integer, the dominant term is exact. 30 * m has at most 10 digits (for subnormal numbers), 31 * and log2_hi has 11 trailing zero bits. 32 * 33 * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h 34 * logF_hi[] + 512 is exact. 35 * 36 * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ... 37 * the leading term is calculated to extra precision in two 38 * parts, the larger of which adds exactly to the dominant 39 * m and F terms. 40 * There are two cases: 41 * 1. when m, j are non-zero (m | j), use absolute 42 * precision for the leading term. 43 * 2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1). 44 * In this case, use a relative precision of 24 bits. 45 * (This is done differently in the original paper) 46 * 47 * Special cases: 48 * 0 return signalling -Inf 49 * neg return signalling NaN 50 * +Inf return +Inf 51 */ 52 53 #if defined(vax) || defined(tahoe) 54 #define _IEEE 0 55 #define TRUNC(x) x = (double) (float) (x) 56 #else 57 #define _IEEE 1 58 #define endian (((*(int *) &one)) ? 1 : 0) 59 #define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000 60 #define infnan(x) 0.0 61 #endif 62 63 double 64 #ifdef _ANSI_SOURCE 65 log(double x) 66 #else 67 log(x) double x; 68 #endif 69 { 70 int m, j; 71 double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0; 72 double logb(), ldexp(); 73 volatile double u1; 74 75 /* Catch special cases */ 76 if (x <= 0) 77 if (_IEEE && x == zero) /* log(0) = -Inf */ 78 return (-one/zero); 79 else if (_IEEE) /* log(neg) = NaN */ 80 return (zero/zero); 81 else if (x == zero) /* NOT REACHED IF _IEEE */ 82 return (infnan(-ERANGE)); 83 else 84 return (infnan(EDOM)); 85 else if (!finite(x)) 86 if (_IEEE) /* x = NaN, Inf */ 87 return (x+x); 88 else 89 return (infnan(ERANGE)); 90 91 /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 92 /* y = F*(1 + f/F) for |f| <= 2^-8 */ 93 94 m = logb(x); 95 g = ldexp(x, -m); 96 if (_IEEE && m == -1022) { 97 j = logb(g), m += j; 98 g = ldexp(g, -j); 99 } 100 j = N*(g-1) + .5; 101 F = (1.0/N) * j + 1; /* F*128 is an integer in [128, 512] */ 102 f = g - F; 103 104 /* Approximate expansion for log(1+f/F) ~= u + q */ 105 g = 1/(2*F+f); 106 u = 2*f*g; 107 v = u*u; 108 q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 109 110 /* case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8, 111 * u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits. 112 * It also adds exactly to |m*log2_hi + log_F_head[j] | < 750 113 */ 114 if (m | j) 115 u1 = u + 513, u1 -= 513; 116 117 /* case 2: |1-x| < 1/256. The m- and j- dependent terms are zero; 118 * u1 = u to 24 bits. 119 */ 120 else 121 u1 = u, TRUNC(u1); 122 u2 = (2.0*(f - F*u1) - u1*f) * g; 123 /* u1 + u2 = 2f/(2F+f) to extra precision. */ 124 125 /* log(x) = log(2^m*F*(1+f/F)) = */ 126 /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); */ 127 /* (exact) + (tiny) */ 128 129 u1 += m*logF_head[N] + logF_head[j]; /* exact */ 130 u2 = (u2 + logF_tail[j]) + q; /* tiny */ 131 u2 += logF_tail[N]*m; 132 return (u1 + u2); 133 } 134 135 /* Extra precision variant, returning 136 * struct {double a, b;}; log(x) = a+b to 63 bits, with 137 * a is rounded to 26 bits. 138 */ 139 struct Double 140 #ifdef _ANSI_SOURCE 141 log__D(double x) 142 #else 143 log__D(x) double x; 144 #endif 145 { 146 int m, j; 147 double F, f, g, q, u, v, u2; 148 double logb(), ldexp(); 149 volatile double u1; 150 struct Double r; 151 152 /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 153 /* y = F*(1 + f/F) for |f| <= 2^-8 */ 154 155 m = logb(x); 156 g = ldexp(x, -m); 157 if (_IEEE && m == -1022) { 158 j = logb(g), m += j; 159 g = ldexp(g, -j); 160 } 161 j = N*(g-1) + .5; 162 F = (1.0/N) * j + 1; 163 f = g - F; 164 165 g = 1/(2*F+f); 166 u = 2*f*g; 167 v = u*u; 168 q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 169 if (m | j) 170 u1 = u + 513, u1 -= 513; 171 else 172 u1 = u, TRUNC(u1); 173 u2 = (2.0*(f - F*u1) - u1*f) * g; 174 175 u1 += m*logF_head[N] + logF_head[j]; 176 177 u2 += logF_tail[j]; u2 += q; 178 u2 += logF_tail[N]*m; 179 r.a = u1 + u2; /* Only difference is here */ 180 TRUNC(r.a); 181 r.b = (u1 - r.a) + u2; 182 return (r); 183 } 184