Refactor quo, exp; improve ln, log10, pow

Various functions (quo, ln, log10, pow, exp) have poor performance
because they made some arbitrary increases to precision. This commit
removes many of those increases and makes some refactors to produce
more accurate results faster.

Quo is now based on the description given on the GDA site. In addition
to being faster, it is not much more accurate, since for some very
bad inputs it would return 0 instead of a correct result. This was
because the c.Precision*2 + 8 precision calculation would be not
precision enough.

Exp is now based on the algorithm given in Variable Precision
Exponential Function by T. E. Hull and A. Abrham (ACM Transactions
on Mathematical Software, Vol 12 #2, pp79–91, ACM, June 1986).

The other functions were able to remove some of their various precision
increases (but not all, see Ln) with great speed increases.

context.go

@@ -15,6 +15,7 @@
 package apd
 
 import (
+	"math"
 	"math/big"
 
 	"github.com/pkg/errors"
@@ -45,6 +46,8 @@ const (
 	// DefaultTraps is the default trap set used by BaseContext. It traps all
 	// flags except Inexact and Rounded.
 	DefaultTraps = ^Condition(0) & ^defaultIgnoreTraps
+
+	errZeroPrecision = "Context may not have 0 Precision for this operation"
 )
 
 // BaseContext is a useful default Context.
@@ -120,7 +123,10 @@ func (c *Context) Quo(d, x, y *Decimal) (Condition, error) {
 	if c.Precision == 0 {
 		// 0 precision is disallowed because we compute the required number of digits
 		// during the 10**x calculation using the precision.
-		return 0, errors.New("Quo requires a Context with > 0 Precision")
+		return 0, errors.New(errZeroPrecision)
+	}
+	if c.Precision > 5000 {
+		return 0, errors.New("Quo requires Precision <= 5000")
 	}
 
 	if y.Coeff.Sign() == 0 {
@@ -133,23 +139,80 @@ func (c *Context) Quo(d, x, y *Decimal) (Condition, error) {
 		}
 		return res.GoError(c.Traps)
 	}
-	a, b, _, err := upscale(x, y)
-	if err != nil {
-		return 0, errors.Wrap(err, "Quo")
+	// An integer variable, adjust, is initialized to 0.
+	var adjust int64
+	// The result coefficient is initialized to 0.
+	quo := new(Decimal)
+	if x.Coeff.Sign() != 0 {
+		dividend := new(big.Int).Abs(&x.Coeff)
+		divisor := new(big.Int).Abs(&y.Coeff)
+
+		// The operand coefficients are adjusted so that the coefficient of the
+		// dividend is greater than or equal to the coefficient of the divisor and
+		// is also less than ten times the coefficient of the divisor, thus:
+
+		// While the coefficient of the dividend is less than the coefficient of
+		// the divisor it is multiplied by 10 and adjust is incremented by 1.
+		for dividend.Cmp(divisor) < 0 {
+			dividend.Mul(dividend, bigTen)
+			adjust++
+		}
+
+		// While the coefficient of the dividend is greater than or equal to ten
+		// times the coefficient of the divisor the coefficient of the divisor is
+		// multiplied by 10 and adjust is decremented by 1.
+		for tmp := new(big.Int); ; {
+			tmp.Mul(divisor, bigTen)
+			if dividend.Cmp(tmp) < 0 {
+				break
+			}
+			divisor.Set(tmp)
+			adjust--
+		}
+
+		// Stage 4, as listed on the GDA site, is to save the dividend to use for
+		// rounding. Instead, we add 3 digits to the desired precision which rounding
+		// will remove. 3 was chosen because it allows all of the non-extended tests
+		// to pass, however there are probably some sets of inputs or contexts that
+		// will produce results with 1ulp of error.
+		prec := int64(c.Precision) + 3
+
+		// The following steps are then repeated until the division is complete:
+		for {
+			// While the coefficient of the divisor is smaller than or equal to the
+			// coefficient of the dividend the former is subtracted from the latter and
+			// the coefficient of the result is incremented by 1.
+			for divisor.Cmp(dividend) <= 0 {
+				dividend.Sub(dividend, divisor)
+				quo.Coeff.Add(&quo.Coeff, bigOne)
+			}
+
+			// If the coefficient of the dividend is now 0 and adjust is greater than
+			// or equal to 0, or if the coefficient of the result has precision digits,
+			// the division is complete.
+			// Step 4, with modifications to the precision as documented above.
+			if (dividend.Sign() == 0 && adjust >= 0) || quo.NumDigits() == prec {
+				break
+			}
+
+			// Otherwise, the coefficients of the result and the dividend are multiplied
+			// by 10 and adjust is incremented by 1.
+			quo.Coeff.Mul(&quo.Coeff, bigTen)
+			dividend.Mul(dividend, bigTen)
+			adjust++
+		}
 	}
 
-	// In order to compute the decimal remainder part, add enough 0s to the
-	// numerator to accurately round with the given precision.
-	// TODO(mjibson): determine a better algorithm for this instead of p*2+8.
-	nc := BaseContext.WithPrecision(c.Precision*2 + 8)
-	f, e, err := exp10(int64(nc.Precision))
-	if err != nil {
-		return 0, nil
+	// The exponent of the result is computed by subtracting the sum of the
+	// original exponent of the divisor and the value of adjust at the end of
+	// the coefficient calculation from the original exponent of the dividend.
+	res := quo.setExponent(c, int64(x.Exponent), int64(-y.Exponent), -adjust)
+
+	// The sign of the result is the exclusive or of the signs of the operands.
+	if xn, yn := x.Sign() == -1, y.Sign() == -1; xn != yn {
+		quo.Coeff.Neg(&quo.Coeff)
 	}
-	f.Mul(a, e)
-	d.Coeff.Quo(f, b)
-	res := d.setExponent(c, -int64(nc.Precision))
-	res |= c.round(d, d)
+	res |= c.round(d, quo)
 	return res.GoError(c.Traps)
 }
 
@@ -381,7 +444,6 @@ func (c *Context) Ln(d, x *Decimal) (Condition, error) {
 	// of how many square-rootings were done using fact and multiply at the end.
 	xr.Set(x)
 	for xr.Cmp(decimalZeroPtNine) < 0 || xr.Cmp(decimalOnePtOne) > 0 {
-		nc.Precision += p
 		ed.Sqrt(xr, xr)
 		ed.Mul(fact, fact, decimalTwo)
 	}
@@ -449,10 +511,15 @@ func (c *Context) Log10(d, x *Decimal) (Condition, error) {
 		return res.GoError(c.Traps)
 	}
 
+	if x.Cmp(decimalOne) == 0 {
+		d.Set(decimalZero)
+		return 0, nil
+	}
+
 	// TODO(mjibson): This is exact under some conditions.
 	res := Inexact
 
-	nc := BaseContext.WithPrecision(c.Precision * 2)
+	nc := BaseContext.WithPrecision(c.Precision + 2)
 	z := new(Decimal)
 	_, err := nc.Ln(z, x)
 	if err != nil {
@@ -467,80 +534,96 @@ func (c *Context) Log10(d, x *Decimal) (Condition, error) {
 	return res.GoError(c.Traps)
 }
 
-// Exp sets d = e**n.
-func (c *Context) Exp(d, n *Decimal) (Condition, error) {
-	if n.Coeff.Sign() == 0 {
+// Exp sets d = e**x.
+func (c *Context) Exp(d, x *Decimal) (Condition, error) {
+	// See: http://dl.acm.org/citation.cfm?id=6498
+
+	if x.Coeff.Sign() == 0 {
 		d.Set(decimalOne)
 		return 0, nil
 	}
 
-	// We are computing (e^n) by splitting n into an integer and a float (e.g
-	// 3.1 ==> x = 3, y = 0.1), this allows us to write e^n = e^(x+y) = e^x * e^y
+	if c.Precision == 0 {
+		return 0, errors.New(errZeroPrecision)
+	}
 
-	integ := new(Decimal)
-	frac := new(Decimal)
-	n.Modf(integ, frac)
+	nc := c.WithPrecision(c.Precision)
+	nc.Rounding = RoundHalfEven
+	res := Inexact | Rounded
 
-	if integ.Exponent > 0 {
-		_, e, err := exp10(int64(integ.Exponent))
-		if err != nil {
-			return 0, err
+	// Stage 1
+	cp := int64(c.Precision)
+	tmp1 := new(Decimal).Abs(x)
+	tmp2 := New(cp*23, 0)
+	// if abs(x) > 23*currentprecision; assert false
+	if tmp1.Cmp(tmp2) > 0 {
+		res |= Overflow
+		if x.Sign() < 0 {
+			res = res.negateOverflowFlags()
 		}
-		integ.Coeff.Mul(&integ.Coeff, e)
-		integ.Exponent = 0
+		return res.GoError(c.Traps)
+	}
+	// if abs(x) <= setexp(.9, -currentprecision); then result 1
+	tmp2.SetCoefficient(9).SetExponent(int32(-cp) - 1)
+	if tmp1.Cmp(tmp2) <= 0 {
+		d.Set(decimalOne)
+		return res.GoError(c.Traps)
 	}
 
-	z := new(Decimal)
-	nc := BaseContext.WithPrecision(c.Precision * 2)
-	nres, err := nc.integerPower(z, decimalE, &integ.Coeff)
-	if err != nil {
-		return nres, errors.Wrap(err, "integerPower")
+	// Stage 2
+	// Add x.NumDigits because the paper assumes that x.Coeff [0.1, 1).
+	t := x.Exponent + int32(x.NumDigits())
+	if t < 0 {
+		t = 0
 	}
-	res := Inexact
-	// TODO(mjibson): figure out a more consistent way to transfer flags from
-	// inner contexts.
-	// These flags can occur here, so transfer that to the parent context.
-	res |= nres & (Overflow | Underflow | Subnormal)
-	return c.smallExp(res, d, z, frac)
-}
+	k := New(1, t)
+	r := new(Decimal)
+	if _, err := nc.Quo(r, x, k); err != nil {
+		return 0, errors.Wrap(err, "QuoInteger")
+	}
+	ra := new(Decimal).Abs(r)
+	p := cp + int64(t) + 2
 
-// smallExp sets d = x * e**y. It should be used with small y values only
-// (|y| < 1).
-func (c *Context) smallExp(res Condition, d, x, y *Decimal) (Condition, error) {
-	n := new(Decimal)
-	e := x.Exponent
-	if e < 0 {
-		e = -e
+	// Stage 3
+	rf, err := ra.Float64()
+	if err != nil {
+		return 0, errors.Wrap(err, "r.Float64")
 	}
-	nc := BaseContext.WithPrecision((uint32(x.NumDigits()) + uint32(e)))
-	if p := c.Precision * 2; nc.Precision < p {
-		nc.Precision = p
+	pf := float64(p)
+	nf := math.Ceil((1.435*pf - 1.182) / math.Log10(pf/rf))
+	if nf > 1000 || math.IsNaN(nf) {
+		return 0, errors.New("too many iterations")
 	}
+	n := int64(nf)
+
+	// Stage 4
+	nc.Precision = uint32(p)
 	ed := NewErrDecimal(nc)
-	z := d
-	tmp := new(Decimal)
-	z.Set(x)
-	tmp.Set(x)
-	for loop := nc.newLoop("exp", z, 1); ; {
-		if err := ed.Err(); err != nil {
-			break
-		}
-		if done, err := loop.done(z); err != nil {
-			return 0, err
-		} else if done {
-			break
-		}
-		ed.Add(n, n, decimalOne)
-		ed.Mul(tmp, tmp, y)
-		ed.Quo(tmp, tmp, n)
-		ed.Add(z, z, tmp)
-	}
-	if err := ed.Err(); err != nil {
+	sum := New(1, 0)
+	tmp2.Exponent = 0
+	for i := n - 1; i > 0; i-- {
+		tmp2.SetCoefficient(i)
+		// tmp1 = r / i
+		ed.Quo(tmp1, r, tmp2)
+		// sum = sum * r / i
+		ed.Mul(sum, tmp1, sum)
+		// sum = sum + 1
+		ed.Add(sum, sum, decimalOne)
+	}
+	if err != ed.Err() {
 		return 0, err
 	}
 
-	// TODO(mjibson): force RoundHalfEven here.
-	res |= c.round(d, z)
+	// sum ** k
+	_, ki, err := exp10(int64(t))
+	if err != nil {
+		return 0, errors.Wrap(err, "ki")
+	}
+	if _, err := nc.integerPower(d, sum, ki); err != nil {
+		return 0, errors.Wrap(err, "integer power")
+	}
+	nc.Precision = c.Precision
+	res |= nc.round(d, d)
 	return res.GoError(c.Traps)
 }
 

decimal.go

@@ -431,3 +431,10 @@ func (d *Decimal) Neg(x *Decimal) *Decimal {
 	d.Coeff.Neg(&d.Coeff)
 	return d
 }
+
+// Abs sets d to |x| and returns d.
+func (d *Decimal) Abs(x *Decimal) *Decimal {
+	d.Set(x)
+	d.Coeff.Abs(&d.Coeff)
+	return d
+}

decimal_test.go

@@ -268,7 +268,7 @@ func TestQuoErr(t *testing.T) {
 		p    uint32
 		err  string
 	}{
-		{x: "1", y: "1", p: 0, err: "Quo requires a Context with > 0 Precision"},
+		{x: "1", y: "1", p: 0, err: errZeroPrecision},
 		{x: "1", y: "0", p: 1, err: "division by zero"},
 	}
 	for _, tc := range tests {