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[] = "@(#)log1p.c 5.2 (Berkeley) 04/29/88"; 20 #endif /* not lint */ 21 22 /* LOG1P(x) 23 * RETURN THE LOGARITHM OF 1+x 24 * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS) 25 * CODED IN C BY K.C. NG, 1/19/85; 26 * REVISED BY K.C. NG on 2/6/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 * 1+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(1+x) = k*ln2 + log(1+f). 53 * 54 * Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers 55 * n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last 56 * 20 bits (for VAX D format), or the last 21 bits ( for IEEE 57 * double) is 0. This ensures n*ln2hi is exactly representable. 58 * 2. In step 1, f may not be representable. A correction term c 59 * for f is computed. It follows that the correction term for 60 * f - t (the leading term of log(1+f) in step 2) is c-c*x. We 61 * add this correction term to n*ln2lo to attenuate the error. 62 * 63 * 64 * Special cases: 65 * log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal; 66 * log1p(INF) is +INF; log1p(-1) is -INF with signal; 67 * only log1p(0)=0 is exact for finite argument. 68 * 69 * Accuracy: 70 * log1p(x) returns the exact log(1+x) nearly rounded. In a test run 71 * with 1,536,000 random arguments on a VAX, the maximum observed 72 * error was .846 ulps (units in the last place). 73 * 74 * Constants: 75 * The hexadecimal values are the intended ones for the following constants. 76 * The decimal values may be used, provided that the compiler will convert 77 * from decimal to binary accurately enough to produce the hexadecimal values 78 * shown. 79 */ 80 81 #if defined(vax)||defined(tahoe) /* VAX D format */ 82 #include <errno.h> 83 #ifdef vax 84 #define _0x(A,B) 0x/**/A/**/B 85 #else /* vax */ 86 #define _0x(A,B) 0x/**/B/**/A 87 #endif /* vax */ 88 /* static double */ 89 /* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */ 90 /* ln2lo = 1.6465949582897081279E-12 , Hex 2^-39 * .E7BCD5E4F1D9CC */ 91 /* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */ 92 static long ln2hix[] = { _0x(7217,4031), _0x(0000,f7d0)}; 93 static long ln2lox[] = { _0x(bcd5,2ce7), _0x(d9cc,e4f1)}; 94 static long sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)}; 95 #define ln2hi (*(double*)ln2hix) 96 #define ln2lo (*(double*)ln2lox) 97 #define sqrt2 (*(double*)sqrt2x) 98 #else /* defined(vax)||defined(tahoe) */ 99 static double 100 ln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */ 101 ln2lo = 1.9082149292705877000E-10 , /*Hex 2^-33 * 1.A39EF35793C76 */ 102 sqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */ 103 #endif /* defined(vax)||defined(tahoe) */ 104 105 double log1p(x) 106 double x; 107 { 108 static double zero=0.0, negone= -1.0, one=1.0, 109 half=1.0/2.0, small=1.0E-20; /* 1+small == 1 */ 110 double logb(),copysign(),scalb(),log__L(),z,s,t,c; 111 int k,finite(); 112 113 #if !defined(vax)&&!defined(tahoe) 114 if(x!=x) return(x); /* x is NaN */ 115 #endif /* !defined(vax)&&!defined(tahoe) */ 116 117 if(finite(x)) { 118 if( x > negone ) { 119 120 /* argument reduction */ 121 if(copysign(x,one)<small) return(x); 122 k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k); 123 if(z+t >= sqrt2 ) 124 { k += 1 ; z *= half; t *= half; } 125 t += negone; x = z + t; 126 c = (t-x)+z ; /* correction term for x */ 127 128 /* compute log(1+x) */ 129 s = x/(2+x); t = x*x*half; 130 c += (k*ln2lo-c*x); 131 z = c+s*(t+log__L(s*s)); 132 x += (z - t) ; 133 134 return(k*ln2hi+x); 135 } 136 /* end of if (x > negone) */ 137 138 else { 139 #if defined(vax)||defined(tahoe) 140 extern double infnan(); 141 if ( x == negone ) 142 return (infnan(-ERANGE)); /* -INF */ 143 else 144 return (infnan(EDOM)); /* NaN */ 145 #else /* defined(vax)||defined(tahoe) */ 146 /* x = -1, return -INF with signal */ 147 if ( x == negone ) return( negone/zero ); 148 149 /* negative argument for log, return NaN with signal */ 150 else return ( zero / zero ); 151 #endif /* defined(vax)||defined(tahoe) */ 152 } 153 } 154 /* end of if (finite(x)) */ 155 156 /* log(-INF) is NaN */ 157 else if(x<0) 158 return(zero/zero); 159 160 /* log(+INF) is INF */ 161 else return(x); 162 } 163