1da2e3ebdSchin #include "FEATURE/uwin"
2da2e3ebdSchin
3da2e3ebdSchin #if !_UWIN || _lib_log1p
4da2e3ebdSchin
_STUB_log1p()5da2e3ebdSchin void _STUB_log1p(){}
6da2e3ebdSchin
7da2e3ebdSchin #else
8da2e3ebdSchin
9da2e3ebdSchin /*
10da2e3ebdSchin * Copyright (c) 1985, 1993
11da2e3ebdSchin * The Regents of the University of California. All rights reserved.
12da2e3ebdSchin *
13da2e3ebdSchin * Redistribution and use in source and binary forms, with or without
14da2e3ebdSchin * modification, are permitted provided that the following conditions
15da2e3ebdSchin * are met:
16da2e3ebdSchin * 1. Redistributions of source code must retain the above copyright
17da2e3ebdSchin * notice, this list of conditions and the following disclaimer.
18da2e3ebdSchin * 2. Redistributions in binary form must reproduce the above copyright
19da2e3ebdSchin * notice, this list of conditions and the following disclaimer in the
20da2e3ebdSchin * documentation and/or other materials provided with the distribution.
21da2e3ebdSchin * 3. Neither the name of the University nor the names of its contributors
22da2e3ebdSchin * may be used to endorse or promote products derived from this software
23da2e3ebdSchin * without specific prior written permission.
24da2e3ebdSchin *
25da2e3ebdSchin * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
26da2e3ebdSchin * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27da2e3ebdSchin * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
28da2e3ebdSchin * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
29da2e3ebdSchin * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30da2e3ebdSchin * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
31da2e3ebdSchin * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
32da2e3ebdSchin * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
33da2e3ebdSchin * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
34da2e3ebdSchin * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
35da2e3ebdSchin * SUCH DAMAGE.
36da2e3ebdSchin */
37da2e3ebdSchin
38da2e3ebdSchin #ifndef lint
39da2e3ebdSchin static char sccsid[] = "@(#)log1p.c 8.1 (Berkeley) 6/4/93";
40da2e3ebdSchin #endif /* not lint */
41da2e3ebdSchin
42da2e3ebdSchin /* LOG1P(x)
43da2e3ebdSchin * RETURN THE LOGARITHM OF 1+x
44da2e3ebdSchin * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS)
45da2e3ebdSchin * CODED IN C BY K.C. NG, 1/19/85;
46da2e3ebdSchin * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85.
47da2e3ebdSchin *
48da2e3ebdSchin * Required system supported functions:
49da2e3ebdSchin * scalb(x,n)
50da2e3ebdSchin * copysign(x,y)
51da2e3ebdSchin * logb(x)
52da2e3ebdSchin * finite(x)
53da2e3ebdSchin *
54da2e3ebdSchin * Required kernel function:
55da2e3ebdSchin * log__L(z)
56da2e3ebdSchin *
57da2e3ebdSchin * Method :
58da2e3ebdSchin * 1. Argument Reduction: find k and f such that
59da2e3ebdSchin * 1+x = 2^k * (1+f),
60da2e3ebdSchin * where sqrt(2)/2 < 1+f < sqrt(2) .
61da2e3ebdSchin *
62da2e3ebdSchin * 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
63da2e3ebdSchin * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
64da2e3ebdSchin * log(1+f) is computed by
65da2e3ebdSchin *
66da2e3ebdSchin * log(1+f) = 2s + s*log__L(s*s)
67da2e3ebdSchin * where
68da2e3ebdSchin * log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
69da2e3ebdSchin *
70da2e3ebdSchin * See log__L() for the values of the coefficients.
71da2e3ebdSchin *
72da2e3ebdSchin * 3. Finally, log(1+x) = k*ln2 + log(1+f).
73da2e3ebdSchin *
74da2e3ebdSchin * Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers
75da2e3ebdSchin * n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last
76da2e3ebdSchin * 20 bits (for VAX D format), or the last 21 bits ( for IEEE
77da2e3ebdSchin * double) is 0. This ensures n*ln2hi is exactly representable.
78da2e3ebdSchin * 2. In step 1, f may not be representable. A correction term c
79da2e3ebdSchin * for f is computed. It follows that the correction term for
80da2e3ebdSchin * f - t (the leading term of log(1+f) in step 2) is c-c*x. We
81da2e3ebdSchin * add this correction term to n*ln2lo to attenuate the error.
82da2e3ebdSchin *
83da2e3ebdSchin *
84da2e3ebdSchin * Special cases:
85da2e3ebdSchin * log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal;
86da2e3ebdSchin * log1p(INF) is +INF; log1p(-1) is -INF with signal;
87da2e3ebdSchin * only log1p(0)=0 is exact for finite argument.
88da2e3ebdSchin *
89da2e3ebdSchin * Accuracy:
90da2e3ebdSchin * log1p(x) returns the exact log(1+x) nearly rounded. In a test run
91da2e3ebdSchin * with 1,536,000 random arguments on a VAX, the maximum observed
92da2e3ebdSchin * error was .846 ulps (units in the last place).
93da2e3ebdSchin *
94da2e3ebdSchin * Constants:
95da2e3ebdSchin * The hexadecimal values are the intended ones for the following constants.
96da2e3ebdSchin * The decimal values may be used, provided that the compiler will convert
97da2e3ebdSchin * from decimal to binary accurately enough to produce the hexadecimal values
98da2e3ebdSchin * shown.
99da2e3ebdSchin */
100da2e3ebdSchin
101da2e3ebdSchin #include <errno.h>
102da2e3ebdSchin #include "mathimpl.h"
103da2e3ebdSchin
104da2e3ebdSchin vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000)
105da2e3ebdSchin vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
106da2e3ebdSchin vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65)
107da2e3ebdSchin
108da2e3ebdSchin ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000)
109da2e3ebdSchin ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76)
110da2e3ebdSchin ic(sqrt2, 1.4142135623730951455E0, 0, 1.6A09E667F3BCD)
111da2e3ebdSchin
112da2e3ebdSchin #ifdef vccast
113da2e3ebdSchin #define ln2hi vccast(ln2hi)
114da2e3ebdSchin #define ln2lo vccast(ln2lo)
115da2e3ebdSchin #define sqrt2 vccast(sqrt2)
116da2e3ebdSchin #endif
117da2e3ebdSchin
118da2e3ebdSchin extern double log1p(x)
119da2e3ebdSchin double x;
120da2e3ebdSchin {
121da2e3ebdSchin const static double zero=0.0, negone= -1.0, one=1.0,
122da2e3ebdSchin half=1.0/2.0, small=1.0E-20; /* 1+small == 1 */
123da2e3ebdSchin double z,s,t,c;
124da2e3ebdSchin int k;
125da2e3ebdSchin
126da2e3ebdSchin #if !defined(vax)&&!defined(tahoe)
127da2e3ebdSchin if(x!=x) return(x); /* x is NaN */
128da2e3ebdSchin #endif /* !defined(vax)&&!defined(tahoe) */
129da2e3ebdSchin
130da2e3ebdSchin if(finite(x)) {
131da2e3ebdSchin if( x > negone ) {
132da2e3ebdSchin
133da2e3ebdSchin /* argument reduction */
134da2e3ebdSchin if(copysign(x,one)<small) return(x);
135da2e3ebdSchin k=(int)logb(one+x); z=scalb(x,-k); t=scalb(one,-k);
136da2e3ebdSchin if(z+t >= sqrt2 )
137da2e3ebdSchin { k += 1 ; z *= half; t *= half; }
138da2e3ebdSchin t += negone; x = z + t;
139da2e3ebdSchin c = (t-x)+z ; /* correction term for x */
140da2e3ebdSchin
141da2e3ebdSchin /* compute log(1+x) */
142da2e3ebdSchin s = x/(2+x); t = x*x*half;
143da2e3ebdSchin c += (k*ln2lo-c*x);
144da2e3ebdSchin z = c+s*(t+__log__L(s*s));
145da2e3ebdSchin x += (z - t) ;
146da2e3ebdSchin
147da2e3ebdSchin return(k*ln2hi+x);
148da2e3ebdSchin }
149da2e3ebdSchin /* end of if (x > negone) */
150da2e3ebdSchin
151da2e3ebdSchin else {
152da2e3ebdSchin #if defined(vax)||defined(tahoe)
153da2e3ebdSchin if ( x == negone )
154da2e3ebdSchin return (infnan(-ERANGE)); /* -INF */
155da2e3ebdSchin else
156da2e3ebdSchin return (infnan(EDOM)); /* NaN */
157da2e3ebdSchin #else /* defined(vax)||defined(tahoe) */
158da2e3ebdSchin /* x = -1, return -INF with signal */
159da2e3ebdSchin if ( x == negone ) return( negone/zero );
160da2e3ebdSchin
161da2e3ebdSchin /* negative argument for log, return NaN with signal */
162da2e3ebdSchin else return ( zero / zero );
163da2e3ebdSchin #endif /* defined(vax)||defined(tahoe) */
164da2e3ebdSchin }
165da2e3ebdSchin }
166da2e3ebdSchin /* end of if (finite(x)) */
167da2e3ebdSchin
168da2e3ebdSchin /* log(-INF) is NaN */
169da2e3ebdSchin else if(x<0)
170da2e3ebdSchin return(zero/zero);
171da2e3ebdSchin
172da2e3ebdSchin /* log(+INF) is INF */
173da2e3ebdSchin else return(x);
174da2e3ebdSchin }
175da2e3ebdSchin
176da2e3ebdSchin #endif
177