Fix exponential()

value/const.go

@@ -81,6 +81,11 @@ func newF(conf *config.Config) *big.Float {
 	return new(big.Float).SetPrec(conf.FloatPrec())
 }
 
+// newFxP returns a new big.Float with extended precision = conf.FloatPrec() + extra.
+func newFxP(conf *config.Config, extra uint) *big.Float {
+	return new(big.Float).SetPrec(conf.FloatPrec() + extra)
+}
+
 func newFloat(c Context) *big.Float {
 	return newF(c.Config())
 }

value/power.go

@@ -5,6 +5,7 @@
 package value
 
 import (
+	"math"
 	"math/big"
 
 	"robpike.io/ivy/config"
@@ -97,32 +98,96 @@ func floatPower(c Context, bx, bexp BigFloat) Value {
 	return BigFloat{z}
 }
 
-// exponential computes exp(x) using the Taylor series. It converges quickly
-// since we call it with only small values of x.
+// exponential computes exp(x) using the Taylor series exp(x) = ∑(x^n/n!) for n ≥ 0.
 func exponential(conf *config.Config, x *big.Float) *big.Float {
-	// The Taylor series for e**x, exp(x), is 1 + x + x²/2! + x³/3! ...
+	// exp(x) is finite if 0.5 × 2^big.MinExp ≤ exp(x) < 1 × 2^big.MaxExp
+	//   ⇒ log(2) × (big.MinExp-1) ≤ x < log(2) × big.MaxExp
+	// While this function properly handles values of x outside of this range,
+	// exit early on extreme values to prevent long running times and simplify the
+	// bounds check to x.exp-1 < log2(big.MaxExp)
+	exp := x.MantExp(nil)
+	if x.IsInf() || exp > 31 {
+		if x.Sign() < 0 {
+			return floatZero
+		}
+		Errorf("exponential overflow")
+	}
+
+	// The following is based on R. P. Brent, P. Zimmermann, Modern Computer
+	// Arithmetic, Cambridge Monographs on Computational and Applied Mathematics
+	// (No. 18), Cambridge University Press
+	// https://members.loria.fr/PZimmermann/mca/pub226.html
+	//
+	// Argument reduction: bring x in the range [0.5, 1)×2^-k for faster
+	// convergence. This also brings extreme values of x for which exp(x) is
+	// 0 or +Inf into a computable range (i.e. for z=x^-k, ∑(z^n/n!) is finite).
+
+	var invert bool
+	z := newF(conf).Set(x)
+	// For z < 0, compute exp(-z) = 1/exp(z).
+	// This is to prevent alternating signs in the power series terms and avoid
+	// cancellation in the summation, as well as keeping the summation in a
+	// known range after argument reduction (1 <= ∑(z^n/n!) < 1+2^(-k+1)).
+	if z.Signbit() {
+		invert = true
+		z.Neg(z)
+	}
+	// §4.3.1 & §4.4.2 (k ≥ 1)
+	k := int(math.Ceil(math.Sqrt(float64(conf.FloatPrec()))))
+	// Extra precision (§4.4)
+	extra := uint(math.Log(float64(conf.FloatPrec()))) + 1
+	if -k < exp {
+		// -k <= -1 < exp
+		exp += k
+		// 0 ≤ k-1 < exp (condition needed to undo argument reduction)
+		z.SetMantExp(z, -exp)
+		// 2 bits of added precision per multiplication when undoing argument reduction.
+		extra += 2 * uint(exp)
+	}
 
-	xN := newF(conf).Set(x)
-	term := newF(conf)
-	n := newF(conf)
-	nFactorial := newF(conf).SetUint64(1)
-	z := newF(conf).SetInt64(1)
+	n := new(big.Float)
+	t0 := newFxP(conf, extra)
+	t1 := newFxP(conf, extra)
+	term := newFxP(conf, extra).SetUint64(1)
+	sum := newFxP(conf, extra).SetUint64(1)
 
-	for loop := newLoop(conf, "exponential", x, 10); ; { // Big exponentials converge slowly.
-		term.Set(xN)
-		term.Quo(term, nFactorial)
-		z.Add(z, term)
+	// TODO: cannot use loop here since it does not handle the extended precision.
+	// term(n) = term(n-1) × x/n is faster than term(n) = x^n / n! (saves one .Mul)
+	for i := uint64(1); ; i++ {
+		t0.Quo(z, n.SetUint64(i))
+		// term.Mul(term, t) and sum.Add(sum, term) require a temp Float for the
+		// result. Manage that ourselves by using our own temps t0, t1, then swap the
+		// pointers.
+		t1.Mul(term, t0)
+		t1, term = term, t1
+		t1.Add(sum, term)
+		t1, sum = sum, t1
 
-		if loop.done(z) {
+		// If term < 1 ulp, we are done. This check is done after the summation since
+		// sum may still change if term ≥ 0.5 ulp, depending on rounding mode.
+		// term < 1 ulp of sum         ⇒ term < 0.5 × 2^(sum.exp-sum.prec+1)
+		// 0 ≤ term < 1 × 2^term.exp   ⇒ 2^term.exp ≤ 2^(sum.exp-sum.prec)
+		// Because of argument reduction, 1 ≤ sum < 1+2^(-k+1) ⇒ sum.exp == 1
+		if term.Sign() == 0 || term.MantExp(nil) <= sum.MantExp(nil) /* ==1 */ -int(sum.Prec()) {
 			break
 		}
-		// Advance x**index (multiply by x).
-		xN.Mul(xN, x)
-		// Advance n, n!.
-		nFactorial.Mul(nFactorial, n.SetUint64(loop.i+1))
 	}
 
-	return z
+	// Undo argument reduction if exp > 0
+	for range exp {
+		// Prevent temp allocations using the same trick as above
+		t0.Mul(sum, sum)
+		t0, sum = sum, t0
+	}
+
+	if invert {
+		// If sum.IsInf the result will be 0 as intended.
+		return z.Quo(n.SetUint64(1), sum)
+	}
+	if sum.IsInf() {
+		Errorf("exponential overflow")
+	}
+	return z.Set(sum)
 }
 
 // integerPower returns x**exp where exp is an int64 of size <= intBits.