1 /* 2 * Copyright (c) 1985 Regents of the University of California. 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms are permitted 6 * provided that this notice is preserved and that due credit is given 7 * to the University of California at Berkeley. The name of the University 8 * may not be used to endorse or promote products derived from this 9 * software without specific prior written permission. This software 10 * is provided ``as is'' without express or implied warranty. 11 * 12 * All recipients should regard themselves as participants in an ongoing 13 * research project and hence should feel obligated to report their 14 * experiences (good or bad) with these elementary function codes, using 15 * the sendbug(8) program, to the authors. 16 */ 17 18 #ifndef lint 19 static char sccsid[] = "@(#)log.c 5.2 (Berkeley) 04/29/88"; 20 #endif /* not lint */ 21 22 /* LOG(X) 23 * RETURN THE LOGARITHM OF x 24 * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS) 25 * CODED IN C BY K.C. NG, 1/19/85; 26 * REVISED BY K.C. NG on 2/7/85, 3/7/85, 3/24/85, 4/16/85. 27 * 28 * Required system supported functions: 29 * scalb(x,n) 30 * copysign(x,y) 31 * logb(x) 32 * finite(x) 33 * 34 * Required kernel function: 35 * log__L(z) 36 * 37 * Method : 38 * 1. Argument Reduction: find k and f such that 39 * x = 2^k * (1+f), 40 * where sqrt(2)/2 < 1+f < sqrt(2) . 41 * 42 * 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) 43 * = 2s + 2/3 s**3 + 2/5 s**5 + ....., 44 * log(1+f) is computed by 45 * 46 * log(1+f) = 2s + s*log__L(s*s) 47 * where 48 * log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...))) 49 * 50 * See log__L() for the values of the coefficients. 51 * 52 * 3. Finally, log(x) = k*ln2 + log(1+f). (Here n*ln2 will be stored 53 * in two floating point number: n*ln2hi + n*ln2lo, n*ln2hi is exact 54 * since the last 20 bits of ln2hi is 0.) 55 * 56 * Special cases: 57 * log(x) is NaN with signal if x < 0 (including -INF) ; 58 * log(+INF) is +INF; log(0) is -INF with signal; 59 * log(NaN) is that NaN with no signal. 60 * 61 * Accuracy: 62 * log(x) returns the exact log(x) nearly rounded. In a test run with 63 * 1,536,000 random arguments on a VAX, the maximum observed error was 64 * .826 ulps (units in the last place). 65 * 66 * Constants: 67 * The hexadecimal values are the intended ones for the following constants. 68 * The decimal values may be used, provided that the compiler will convert 69 * from decimal to binary accurately enough to produce the hexadecimal values 70 * shown. 71 */ 72 73 #if defined(vax)||defined(tahoe) /* VAX D format */ 74 #include <errno.h> 75 #ifdef vax 76 #define _0x(A,B) 0x/**/A/**/B 77 #else /* vax */ 78 #define _0x(A,B) 0x/**/B/**/A 79 #endif /* vax */ 80 /* static double */ 81 /* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */ 82 /* ln2lo = 1.6465949582897081279E-12 , Hex 2^-39 * .E7BCD5E4F1D9CC */ 83 /* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */ 84 static long ln2hix[] = { _0x(7217,4031), _0x(0000,f7d0)}; 85 static long ln2lox[] = { _0x(bcd5,2ce7), _0x(d9cc,e4f1)}; 86 static long sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)}; 87 #define ln2hi (*(double*)ln2hix) 88 #define ln2lo (*(double*)ln2lox) 89 #define sqrt2 (*(double*)sqrt2x) 90 #else /* defined(vax)||defined(tahoe) */ 91 static double 92 ln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */ 93 ln2lo = 1.9082149292705877000E-10 , /*Hex 2^-33 * 1.A39EF35793C76 */ 94 sqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */ 95 #endif /* defined(vax)||defined(tahoe) */ 96 97 double log(x) 98 double x; 99 { 100 static double zero=0.0, negone= -1.0, half=1.0/2.0; 101 double logb(),scalb(),copysign(),log__L(),s,z,t; 102 int k,n,finite(); 103 104 #if !defined(vax)&&!defined(tahoe) 105 if(x!=x) return(x); /* x is NaN */ 106 #endif /* !defined(vax)&&!defined(tahoe) */ 107 if(finite(x)) { 108 if( x > zero ) { 109 110 /* argument reduction */ 111 k=logb(x); x=scalb(x,-k); 112 if(k == -1022) /* subnormal no. */ 113 {n=logb(x); x=scalb(x,-n); k+=n;} 114 if(x >= sqrt2 ) {k += 1; x *= half;} 115 x += negone ; 116 117 /* compute log(1+x) */ 118 s=x/(2+x); t=x*x*half; 119 z=k*ln2lo+s*(t+log__L(s*s)); 120 x += (z - t) ; 121 122 return(k*ln2hi+x); 123 } 124 /* end of if (x > zero) */ 125 126 else { 127 #if defined(vax)||defined(tahoe) 128 extern double infnan(); 129 if ( x == zero ) 130 return (infnan(-ERANGE)); /* -INF */ 131 else 132 return (infnan(EDOM)); /* NaN */ 133 #else /* defined(vax)||defined(tahoe) */ 134 /* zero argument, return -INF with signal */ 135 if ( x == zero ) 136 return( negone/zero ); 137 138 /* negative argument, return NaN with signal */ 139 else 140 return ( zero / zero ); 141 #endif /* defined(vax)||defined(tahoe) */ 142 } 143 } 144 /* end of if (finite(x)) */ 145 /* NOTREACHED if defined(vax)||defined(tahoe) */ 146 147 /* log(-INF) is NaN with signal */ 148 else if (x<0) 149 return(zero/zero); 150 151 /* log(+INF) is +INF */ 152 else return(x); 153 154 } 155