Branch data Line data Source code
1 : : /*
2 : : * Copyright (c) 1985, 1993
3 : : * The Regents of the University of California. All rights reserved.
4 : : *
5 : : * Redistribution and use in source and binary forms, with or without
6 : : * modification, are permitted provided that the following conditions
7 : : * are met:
8 : : * 1. Redistributions of source code must retain the above copyright
9 : : * notice, this list of conditions and the following disclaimer.
10 : : * 2. Redistributions in binary form must reproduce the above copyright
11 : : * notice, this list of conditions and the following disclaimer in the
12 : : * documentation and/or other materials provided with the distribution.
13 : : * 3. Neither the name of the University nor the names of its contributors
14 : : * may be used to endorse or promote products derived from this software
15 : : * without specific prior written permission.
16 : : *
17 : : * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
18 : : * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 : : * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 : : * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
21 : : * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22 : : * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23 : : * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24 : : * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25 : : * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26 : : * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27 : : * SUCH DAMAGE.
28 : : */
29 : :
30 : : /* @(#)exp.c 8.1 (Berkeley) 6/4/93 */
31 : : #include "cdefs-compat.h"
32 : : //__FBSDID("$FreeBSD: src/lib/msun/bsdsrc/b_exp.c,v 1.9 2011/10/16 05:37:20 das Exp $");
33 : :
34 : : #include <openlibm_math.h>
35 : :
36 : : /* EXP(X)
37 : : * RETURN THE EXPONENTIAL OF X
38 : : * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
39 : : * CODED IN C BY K.C. NG, 1/19/85;
40 : : * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
41 : : *
42 : : * Required system supported functions:
43 : : * scalb(x,n)
44 : : * copysign(x,y)
45 : : * finite(x)
46 : : *
47 : : * Method:
48 : : * 1. Argument Reduction: given the input x, find r and integer k such
49 : : * that
50 : : * x = k*ln2 + r, |r| <= 0.5*ln2 .
51 : : * r will be represented as r := z+c for better accuracy.
52 : : *
53 : : * 2. Compute exp(r) by
54 : : *
55 : : * exp(r) = 1 + r + r*R1/(2-R1),
56 : : * where
57 : : * R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))).
58 : : *
59 : : * 3. exp(x) = 2^k * exp(r) .
60 : : *
61 : : * Special cases:
62 : : * exp(INF) is INF, exp(NaN) is NaN;
63 : : * exp(-INF)= 0;
64 : : * for finite argument, only exp(0)=1 is exact.
65 : : *
66 : : * Accuracy:
67 : : * exp(x) returns the exponential of x nearly rounded. In a test run
68 : : * with 1,156,000 random arguments on a VAX, the maximum observed
69 : : * error was 0.869 ulps (units in the last place).
70 : : */
71 : :
72 : : #include "mathimpl.h"
73 : :
74 : : static const double p1 = 0x1.555555555553ep-3;
75 : : static const double p2 = -0x1.6c16c16bebd93p-9;
76 : : static const double p3 = 0x1.1566aaf25de2cp-14;
77 : : static const double p4 = -0x1.bbd41c5d26bf1p-20;
78 : : static const double p5 = 0x1.6376972bea4d0p-25;
79 : : static const double ln2hi = 0x1.62e42fee00000p-1;
80 : : static const double ln2lo = 0x1.a39ef35793c76p-33;
81 : : static const double lnhuge = 0x1.6602b15b7ecf2p9;
82 : : static const double lntiny = -0x1.77af8ebeae354p9;
83 : : static const double invln2 = 0x1.71547652b82fep0;
84 : :
85 : : #if 0
86 : : OLM_DLLEXPORT double exp(x)
87 : : double x;
88 : : {
89 : : double z,hi,lo,c;
90 : : int k;
91 : :
92 : : #if !defined(vax)&&!defined(tahoe)
93 : : if(x!=x) return(x); /* x is NaN */
94 : : #endif /* !defined(vax)&&!defined(tahoe) */
95 : : if( x <= lnhuge ) {
96 : : if( x >= lntiny ) {
97 : :
98 : : /* argument reduction : x --> x - k*ln2 */
99 : :
100 : : k=invln2*x+copysign(0.5,x); /* k=NINT(x/ln2) */
101 : :
102 : : /* express x-k*ln2 as hi-lo and let x=hi-lo rounded */
103 : :
104 : : hi=x-k*ln2hi;
105 : : x=hi-(lo=k*ln2lo);
106 : :
107 : : /* return 2^k*[1+x+x*c/(2+c)] */
108 : : z=x*x;
109 : : c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
110 : : return scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k);
111 : :
112 : : }
113 : : /* end of x > lntiny */
114 : :
115 : : else
116 : : /* exp(-big#) underflows to zero */
117 : : if(finite(x)) return(scalb(1.0,-5000));
118 : :
119 : : /* exp(-INF) is zero */
120 : : else return(0.0);
121 : : }
122 : : /* end of x < lnhuge */
123 : :
124 : : else
125 : : /* exp(INF) is INF, exp(+big#) overflows to INF */
126 : : return( finite(x) ? scalb(1.0,5000) : x);
127 : : }
128 : : #endif
129 : :
130 : : /* returns exp(r = x + c) for |c| < |x| with no overlap. */
131 : :
132 : 0 : double __exp__D(x, c)
133 : : double x, c;
134 : : {
135 : : double z,hi,lo;
136 : : int k;
137 : :
138 [ # # ]: 0 : if (x != x) /* x is NaN */
139 : 0 : return(x);
140 [ # # ]: 0 : if ( x <= lnhuge ) {
141 [ # # ]: 0 : if ( x >= lntiny ) {
142 : :
143 : : /* argument reduction : x --> x - k*ln2 */
144 : 0 : z = invln2*x;
145 : 0 : k = z + copysign(.5, x);
146 : :
147 : : /* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */
148 : :
149 : 0 : hi=(x-k*ln2hi); /* Exact. */
150 : 0 : x= hi - (lo = k*ln2lo-c);
151 : : /* return 2^k*[1+x+x*c/(2+c)] */
152 : 0 : z=x*x;
153 : 0 : c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
154 : 0 : c = (x*c)/(2.0-c);
155 : :
156 : 0 : return scalbn(1.+(hi-(lo - c)), k);
157 : : }
158 : : /* end of x > lntiny */
159 : :
160 : : else
161 : : /* exp(-big#) underflows to zero */
162 [ # # ]: 0 : if(isfinite(x)) return(scalbn(1.0,-5000));
163 : :
164 : : /* exp(-INF) is zero */
165 : 0 : else return(0.0);
166 : : }
167 : : /* end of x < lnhuge */
168 : :
169 : : else
170 : : /* exp(INF) is INF, exp(+big#) overflows to INF */
171 [ # # ]: 0 : return( isfinite(x) ? scalbn(1.0,5000) : x);
172 : : }
|
