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.7 (Berkeley) 12/02/92"; 10 #endif /* not lint */ 11 12 #include <math.h> 13 #include <errno.h> 14 15 #include "log_table.h" 16 #include "dmath.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) (double) (float) (x) 56 #else 57 #define _IEEE 1 58 #define TRUNC(x) *(((int *) &x) + 1) &= 0xf8000000 59 #define infnan(x) 0.0 60 #endif 61 62 double 63 #ifdef _ANSI_SOURCE 64 log(double x) 65 #else 66 log(x) double x; 67 #endif 68 { 69 int m, j; 70 double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0; 71 double logb(), ldexp(); 72 volatile double u1; 73 74 /* Catch special cases */ 75 if (x <= 0) 76 if (_IEEE && x == zero) /* log(0) = -Inf */ 77 return (-one/zero); 78 else if (_IEEE) /* log(neg) = NaN */ 79 return (zero/zero); 80 else if (x == zero) /* NOT REACHED IF _IEEE */ 81 return (infnan(-ERANGE)); 82 else 83 return (infnan(EDOM)); 84 else if (!finite(x)) 85 if (_IEEE) /* x = NaN, Inf */ 86 return (x+x); 87 else 88 return (infnan(ERANGE)); 89 90 /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 91 /* y = F*(1 + f/F) for |f| <= 2^-8 */ 92 93 m = logb(x); 94 g = ldexp(x, -m); 95 if (_IEEE && m == -1022) { 96 j = logb(g), m += j; 97 g = ldexp(g, -j); 98 } 99 j = N*(g-1) + .5; 100 F = (1.0/N) * j + 1; /* F*128 is an integer in [128, 512] */ 101 f = g - F; 102 103 /* Approximate expansion for log(1+f/F) ~= u + q */ 104 g = 1/(2*F+f); 105 u = 2*f*g; 106 v = u*u; 107 q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 108 109 /* case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8, 110 * u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits. 111 * It also adds exactly to |m*log2_hi + log_F_head[j] | < 750 112 */ 113 if (m | j) 114 u1 = u + 513, u1 -= 513; 115 116 /* case 2: |1-x| < 1/256. The m- and j- dependent terms are zero; 117 * u1 = u to 24 bits. 118 */ 119 else 120 u1 = u, TRUNC(u1); 121 u2 = (2.0*(f - F*u1) - u1*f) * g; 122 /* u1 + u2 = 2f/(2F+f) to extra precision. */ 123 124 /* log(x) = log(2^m*F*(1+f/F)) = */ 125 /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); */ 126 /* (exact) + (tiny) */ 127 128 u1 += m*logF_head[N] + logF_head[j]; /* exact */ 129 u2 = (u2 + logF_tail[j]) + q; /* tiny */ 130 u2 += logF_tail[N]*m; 131 return (u1 + u2); 132 } 133 134 /* Extra precision variant, returning 135 * struct {double a, b;}; log(x) = a+b to 63 bits, with 136 * a is rounded to 26 bits. 137 */ 138 struct Double 139 #ifdef _ANSI_SOURCE 140 log__D(double x) 141 #else 142 log__D(x) double x; 143 #endif 144 { 145 int m, j; 146 double F, f, g, q, u, v, u2; 147 double logb(), ldexp(); 148 volatile double u1; 149 struct Double r; 150 151 /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 152 /* y = F*(1 + f/F) for |f| <= 2^-8 */ 153 154 m = logb(x); 155 g = ldexp(x, -m); 156 if (_IEEE && m == -1022) { 157 j = logb(g), m += j; 158 g = ldexp(g, -j); 159 } 160 j = N*(g-1) + .5; 161 F = (1.0/N) * j + 1; 162 f = g - F; 163 164 g = 1/(2*F+f); 165 u = 2*f*g; 166 v = u*u; 167 q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 168 if (m | j) 169 u1 = u + 513, u1 -= 513; 170 else 171 u1 = u, TRUNC(u1); 172 u2 = (2.0*(f - F*u1) - u1*f) * g; 173 174 u1 += m*logF_head[N] + logF_head[j]; 175 176 u2 += logF_tail[j]; u2 += q; 177 u2 += logF_tail[N]*m; 178 r.a = u1 + u2; /* Only difference is here */ 179 TRUNC(r.a); 180 r.b = (u1 - r.a) + u2; 181 return (r); 182 } 183