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