Source: expm1.go in package math


package math

The original C code, the long comment, and the constants below are from FreeBSD's /usr/src/lib/msun/src/s_expm1.c and came with this notice. The go code is a simplified version of the original C. ==================================================== Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. Developed at SunPro, a Sun Microsystems, Inc. business. Permission to use, copy, modify, and distribute this software is freely granted, provided that this notice is preserved. ==================================================== expm1(x) Returns exp(x)-1, the exponential of x minus 1. Method 1. Argument reduction: Given x, find r and integer k such that x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 Here a correction term c will be computed to compensate the error in r when rounded to a floating-point number. 2. Approximating expm1(r) by a special rational function on the interval [0,0.34658]: Since r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 - r**4/360 + ... we define R1(r*r) by r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 * R1(r*r) That is, R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) = 1 - r**2/60 + r**4/2520 - r**6/100800 + ... We use a special Reme algorithm on [0,0.347] to generate a polynomial of degree 5 in r*r to approximate R1. The maximum error of this polynomial approximation is bounded by 2**-61. In other words, R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 where Q1 = -1.6666666666666567384E-2, Q2 = 3.9682539681370365873E-4, Q3 = -9.9206344733435987357E-6, Q4 = 2.5051361420808517002E-7, Q5 = -6.2843505682382617102E-9; (where z=r*r, and the values of Q1 to Q5 are listed below) with error bounded by | 5 | -61 | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 | | expm1(r) = exp(r)-1 is then computed by the following specific way which minimize the accumulation rounding error: 2 3 r r [ 3 - (R1 + R1*r/2) ] expm1(r) = r + --- + --- * [--------------------] 2 2 [ 6 - r*(3 - R1*r/2) ] To compensate the error in the argument reduction, we use expm1(r+c) = expm1(r) + c + expm1(r)*c ~ expm1(r) + c + r*c Thus c+r*c will be added in as the correction terms for expm1(r+c). Now rearrange the term to avoid optimization screw up: ( 2 2 ) ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) ( ) = r - E 3. Scale back to obtain expm1(x): From step 1, we have expm1(x) = either 2**k*[expm1(r)+1] - 1 = or 2**k*[expm1(r) + (1-2**-k)] 4. Implementation notes: (A). To save one multiplication, we scale the coefficient Qi to Qi*2**i, and replace z by (x**2)/2. (B). To achieve maximum accuracy, we compute expm1(x) by (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf) (ii) if k=0, return r-E (iii) if k=-1, return 0.5*(r-E)-0.5 (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) else return 1.0+2.0*(r-E); (v) if (k<-2||k>56) return 2**k(1-(E-r)) - 1 (or exp(x)-1) (vi) if k <= 20, return 2**k((1-2**-k)-(E-r)), else (vii) return 2**k(1-((E+2**-k)-r)) Special cases: expm1(INF) is INF, expm1(NaN) is NaN; expm1(-INF) is -1, and for finite argument, only expm1(0)=0 is exact. Accuracy: according to an error analysis, the error is always less than 1 ulp (unit in the last place). Misc. info. For IEEE double if x > 7.09782712893383973096e+02 then expm1(x) overflow Constants: The hexadecimal values are the intended ones for the following constants. The decimal values may be used, provided that the compiler will convert from decimal to binary accurately enough to produce the hexadecimal values shown.

Expm1 returns e**x - 1, the base-e exponential of x minus 1. It is more accurate than Exp(x) - 1 when x is near zero. Special cases are: Expm1(+Inf) = +Inf Expm1(-Inf) = -1 Expm1(NaN) = NaN Very large values overflow to -1 or +Inf.

func Expm1(x float64) float64

func expm1(x float64) float64 {
	const (
		Othreshold = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF
		Ln2X56     = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1
		Ln2HalfX3  = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73
		Ln2Half    = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef
		Ln2Hi      = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
		Ln2Lo      = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
		InvLn2     = 1.44269504088896338700e+00 // 0x3ff71547652b82fe

scaled coefficients related to expm1

		Q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4
		Q2 = 1.58730158725481460165e-03  // 0x3F5A01A019FE5585
		Q3 = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7
		Q4 = 4.00821782732936239552e-06  // 0x3ED0CFCA86E65239
		Q5 = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D
	)

special cases

	switch {
	case IsInf(x, 1) || IsNaN(x):
		return x
	case IsInf(x, -1):
		return -1
	}

	absx := x
	sign := false
	if x < 0 {
		absx = -absx
		sign = true
	}

filter out huge argument

	if absx >= Ln2X56 { // if |x| >= 56 * ln2
		if sign {
			return -1 // x < -56*ln2, return -1
		}
		if absx >= Othreshold { // if |x| >= 709.78...
			return Inf(1)
		}
	}

argument reduction

	var c float64
	var k int
	if absx > Ln2Half { // if  |x| > 0.5 * ln2
		var hi, lo float64
		if absx < Ln2HalfX3 { // and |x| < 1.5 * ln2
			if !sign {
				hi = x - Ln2Hi
				lo = Ln2Lo
				k = 1
			} else {
				hi = x + Ln2Hi
				lo = -Ln2Lo
				k = -1
			}
		} else {
			if !sign {
				k = int(InvLn2*x + 0.5)
			} else {
				k = int(InvLn2*x - 0.5)
			}
			t := float64(k)
			hi = x - t*Ln2Hi // t * Ln2Hi is exact here
			lo = t * Ln2Lo
		}
		x = hi - lo
		c = (hi - x) - lo
	} else if absx < Tiny { // when |x| < 2**-54, return x
		return x
	} else {
		k = 0
	}

x is now in primary range

	hfx := 0.5 * x
	hxs := x * hfx
	r1 := 1 + hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))))
	t := 3 - r1*hfx
	e := hxs * ((r1 - t) / (6.0 - x*t))
	if k == 0 {
		return x - (x*e - hxs) // c is 0
	}
	e = (x*(e-c) - c)
	e -= hxs
	switch {
	case k == -1:
		return 0.5*(x-e) - 0.5
	case k == 1:
		if x < -0.25 {
			return -2 * (e - (x + 0.5))
		}
		return 1 + 2*(x-e)
	case k <= -2 || k > 56: // suffice to return exp(x)-1
		y := 1 - (e - x)
		y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
		return y - 1
	}
	if k < 20 {
		t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2**-k
		y := t - (e - x)
		y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
		return y
	}
	t = Float64frombits(uint64(0x3ff-k) << 52) // 2**-k
	y := x - (e + t)
	y++
	y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
	return y