1 // Copyright 2010 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
4 5 package math
6 7 // The original C code, the long comment, and the constants
8 // below are from FreeBSD's /usr/src/lib/msun/src/s_expm1.c
9 // and came with this notice. The go code is a simplified
10 // version of the original C.
11 //
12 // ====================================================
13 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
14 //
15 // Developed at SunPro, a Sun Microsystems, Inc. business.
16 // Permission to use, copy, modify, and distribute this
17 // software is freely granted, provided that this notice
18 // is preserved.
19 // ====================================================
20 //
21 // expm1(x)
22 // Returns exp(x)-1, the exponential of x minus 1.
23 //
24 // Method
25 // 1. Argument reduction:
26 // Given x, find r and integer k such that
27 //
28 // x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658
29 //
30 // Here a correction term c will be computed to compensate
31 // the error in r when rounded to a floating-point number.
32 //
33 // 2. Approximating expm1(r) by a special rational function on
34 // the interval [0,0.34658]:
35 // Since
36 // r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 - r**4/360 + ...
37 // we define R1(r*r) by
38 // r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 * R1(r*r)
39 // That is,
40 // R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
41 // = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
42 // = 1 - r**2/60 + r**4/2520 - r**6/100800 + ...
43 // We use a special Reme algorithm on [0,0.347] to generate
44 // a polynomial of degree 5 in r*r to approximate R1. The
45 // maximum error of this polynomial approximation is bounded
46 // by 2**-61. In other words,
47 // R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
48 // where Q1 = -1.6666666666666567384E-2,
49 // Q2 = 3.9682539681370365873E-4,
50 // Q3 = -9.9206344733435987357E-6,
51 // Q4 = 2.5051361420808517002E-7,
52 // Q5 = -6.2843505682382617102E-9;
53 // (where z=r*r, and the values of Q1 to Q5 are listed below)
54 // with error bounded by
55 // | 5 | -61
56 // | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2
57 // | |
58 //
59 // expm1(r) = exp(r)-1 is then computed by the following
60 // specific way which minimize the accumulation rounding error:
61 // 2 3
62 // r r [ 3 - (R1 + R1*r/2) ]
63 // expm1(r) = r + --- + --- * [--------------------]
64 // 2 2 [ 6 - r*(3 - R1*r/2) ]
65 //
66 // To compensate the error in the argument reduction, we use
67 // expm1(r+c) = expm1(r) + c + expm1(r)*c
68 // ~ expm1(r) + c + r*c
69 // Thus c+r*c will be added in as the correction terms for
70 // expm1(r+c). Now rearrange the term to avoid optimization
71 // screw up:
72 // ( 2 2 )
73 // ({ ( r [ R1 - (3 - R1*r/2) ] ) } r )
74 // expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
75 // ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 )
76 // ( )
77 //
78 // = r - E
79 // 3. Scale back to obtain expm1(x):
80 // From step 1, we have
81 // expm1(x) = either 2**k*[expm1(r)+1] - 1
82 // = or 2**k*[expm1(r) + (1-2**-k)]
83 // 4. Implementation notes:
84 // (A). To save one multiplication, we scale the coefficient Qi
85 // to Qi*2**i, and replace z by (x**2)/2.
86 // (B). To achieve maximum accuracy, we compute expm1(x) by
87 // (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
88 // (ii) if k=0, return r-E
89 // (iii) if k=-1, return 0.5*(r-E)-0.5
90 // (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E)
91 // else return 1.0+2.0*(r-E);
92 // (v) if (k<-2||k>56) return 2**k(1-(E-r)) - 1 (or exp(x)-1)
93 // (vi) if k <= 20, return 2**k((1-2**-k)-(E-r)), else
94 // (vii) return 2**k(1-((E+2**-k)-r))
95 //
96 // Special cases:
97 // expm1(INF) is INF, expm1(NaN) is NaN;
98 // expm1(-INF) is -1, and
99 // for finite argument, only expm1(0)=0 is exact.
100 //
101 // Accuracy:
102 // according to an error analysis, the error is always less than
103 // 1 ulp (unit in the last place).
104 //
105 // Misc. info.
106 // For IEEE double
107 // if x > 7.09782712893383973096e+02 then expm1(x) overflow
108 //
109 // Constants:
110 // The hexadecimal values are the intended ones for the following
111 // constants. The decimal values may be used, provided that the
112 // compiler will convert from decimal to binary accurately enough
113 // to produce the hexadecimal values shown.
114 //
115 116 // Expm1 returns e**x - 1, the base-e exponential of x minus 1.
117 // It is more accurate than [Exp](x) - 1 when x is near zero.
118 //
119 // Special cases are:
120 //
121 // Expm1(+Inf) = +Inf
122 // Expm1(-Inf) = -1
123 // Expm1(NaN) = NaN
124 //
125 // Very large values overflow to -1 or +Inf.
126 func Expm1(x float64) float64 {
127 if haveArchExpm1 {
128 return archExpm1(x)
129 }
130 return expm1(x)
131 }
132 133 func expm1(x float64) float64 {
134 const (
135 Othreshold = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF
136 Ln2X56 = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1
137 Ln2HalfX3 = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73
138 Ln2Half = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef
139 Ln2Hi = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
140 Ln2Lo = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
141 InvLn2 = 1.44269504088896338700e+00 // 0x3ff71547652b82fe
142 Tiny = 1.0 / (1 << 54) // 2**-54 = 0x3c90000000000000
143 // scaled coefficients related to expm1
144 Q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4
145 Q2 = 1.58730158725481460165e-03 // 0x3F5A01A019FE5585
146 Q3 = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7
147 Q4 = 4.00821782732936239552e-06 // 0x3ED0CFCA86E65239
148 Q5 = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D
149 )
150 151 // special cases
152 switch {
153 case IsInf(x, 1) || IsNaN(x):
154 return x
155 case IsInf(x, -1):
156 return -1
157 }
158 159 absx := x
160 sign := false
161 if x < 0 {
162 absx = -absx
163 sign = true
164 }
165 166 // filter out huge argument
167 if absx >= Ln2X56 { // if |x| >= 56 * ln2
168 if sign {
169 return -1 // x < -56*ln2, return -1
170 }
171 if absx >= Othreshold { // if |x| >= 709.78...
172 return Inf(1)
173 }
174 }
175 176 // argument reduction
177 var c float64
178 var k int
179 if absx > Ln2Half { // if |x| > 0.5 * ln2
180 var hi, lo float64
181 if absx < Ln2HalfX3 { // and |x| < 1.5 * ln2
182 if !sign {
183 hi = x - Ln2Hi
184 lo = Ln2Lo
185 k = 1
186 } else {
187 hi = x + Ln2Hi
188 lo = -Ln2Lo
189 k = -1
190 }
191 } else {
192 if !sign {
193 k = int(InvLn2*x + 0.5)
194 } else {
195 k = int(InvLn2*x - 0.5)
196 }
197 t := float64(k)
198 hi = x - t*Ln2Hi // t * Ln2Hi is exact here
199 lo = t * Ln2Lo
200 }
201 x = hi - lo
202 c = (hi - x) - lo
203 } else if absx < Tiny { // when |x| < 2**-54, return x
204 return x
205 } else {
206 k = 0
207 }
208 209 // x is now in primary range
210 hfx := 0.5 * x
211 hxs := x * hfx
212 r1 := 1 + hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))))
213 t := 3 - r1*hfx
214 e := hxs * ((r1 - t) / (6.0 - x*t))
215 if k == 0 {
216 return x - (x*e - hxs) // c is 0
217 }
218 e = (x*(e-c) - c)
219 e -= hxs
220 switch {
221 case k == -1:
222 return 0.5*(x-e) - 0.5
223 case k == 1:
224 if x < -0.25 {
225 return -2 * (e - (x + 0.5))
226 }
227 return 1 + 2*(x-e)
228 case k <= -2 || k > 56: // suffice to return exp(x)-1
229 y := 1 - (e - x)
230 y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
231 return y - 1
232 }
233 if k < 20 {
234 t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2**-k
235 y := t - (e - x)
236 y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
237 return y
238 }
239 t = Float64frombits(uint64(0x3ff-k) << 52) // 2**-k
240 y := x - (e + t)
241 y++
242 y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
243 return y
244 }
245