diff --git a/context.go b/context.go index 0d3125d..0b804cd 100644 --- a/context.go +++ b/context.go @@ -284,100 +284,85 @@ func (c *Context) Quo(d, x, y *Decimal) (Condition, error) { return res, err } - if c.Precision > 5000 { - // High precision could result in a large number of iterations. Arbitrarily - // limit the precision to prevent runaway processes. This limit was chosen - // arbitrarily and could likely be increased or removed if the impact was - // measured. Until then, this is an attempt to prevent users from shooting - // themselves in the foot. - return 0, errors.New("Quo requires Precision <= 5000") - } - // The sign of the result is the exclusive or of the signs of the operands. neg := x.Negative != y.Negative - // An integer variable, adjust, is initialized to 0. - var adjust int64 - // The result coefficient is initialized to 0. - var quo Decimal - var res Condition - var diff int64 - if !x.IsZero() { - var dividend, divisor BigInt - dividend.Abs(&x.Coeff) - divisor.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(÷nd, bigTen) - adjust++ - } + // Shift the resulting exponent by the difference between the dividend and + // the divisor's exponent after performing arithmetic on the coefficients. + shift := int64(x.Exponent - y.Exponent) - // 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. - var tmp BigInt - for { - tmp.Mul(&divisor, bigTen) - if dividend.Cmp(&tmp) < 0 { - break - } - divisor.Set(&tmp) - adjust-- - } - - prec := int64(c.Precision) - - // 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(÷nd) <= 0 { - dividend.Sub(÷nd, &divisor) - quo.Coeff.Add(&quo.Coeff, bigOne) - } + var res Condition + if x.IsZero() { + d.Set(decimalZero) + d.Negative = neg + res |= d.setExponent(c, unknownNumDigits, res, shift) + return c.goError(res) + } - // 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. - if (dividend.Sign() == 0 && adjust >= 0) || quo.NumDigits() == prec { - break - } + var dividend, divisor BigInt + dividend.Abs(&x.Coeff) + divisor.Abs(&y.Coeff) - // 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(÷nd, bigTen) - adjust++ - } + // 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. While doing + // so, keep track of how far the two have been adjusted. + ndDividend := NumDigits(÷nd) + ndDivisor := NumDigits(&divisor) + ndDiff := ndDividend - ndDivisor + var tmpE BigInt + if ndDiff < 0 { + // numDigits(dividend) < numDigits(divisor), multiply dividend by 10^diff. + dividend.Mul(÷nd, tableExp10(-ndDiff, &tmpE)) + } else if ndDiff > 0 { + // numDigits(dividend) > numDigits(divisor), multiply divisor by 10^diff. + divisor.Mul(&divisor, tableExp10(ndDiff, &tmpE)) + } + adjCoeffs := -ndDiff + if dividend.Cmp(&divisor) < 0 { + // dividend < divisor, multiply dividend by 10. + dividend.Mul(÷nd, bigTen) + adjCoeffs++ + } + + // In order to compute the decimal remainder part, add enough 0s to the + // numerator to accurately round with the given precision. -1 because the + // previous adjustment ensured that the dividend is already greater than or + // equal to the divisor, so the result will always be greater than or equal + // to 1. + adjExp10 := int64(c.Precision - 1) + dividend.Mul(÷nd, tableExp10(adjExp10, &tmpE)) + + // Perform the division. + var rem BigInt + d.Coeff.QuoRem(÷nd, &divisor, &rem) + d.Form = Finite + d.Negative = neg - // Use the adjusted exponent to determine if we are Subnormal. If so, - // don't round. - adj := int64(x.Exponent) + int64(-y.Exponent) - adjust + quo.NumDigits() - 1 - // Any remainder (the final coefficient of the dividend) is recorded and - // taken into account for rounding. - if dividend.Sign() != 0 && adj >= int64(c.MinExponent) { + // If there was a remainder, it is taken into account for rounding. To do + // so, we determine whether the remainder was more or less than half of the + // divisor and round accordingly. + nd := NumDigits(&d.Coeff) + if rem.Sign() != 0 { + // Use the adjusted exponent to determine if we are Subnormal. + // If so, don't round. This computation of adj and the check + // against MinExponent mirrors the logic in setExponent. + adj := shift + (-adjCoeffs) + (-adjExp10) + nd - 1 + if adj >= int64(c.MinExponent) { res |= Inexact | Rounded - dividend.Mul(÷nd, bigTwo) - half := dividend.Cmp(&divisor) - if c.Rounding.ShouldAddOne(&quo.Coeff, quo.Negative, half) { - roundAddOne(&quo.Coeff, &diff) + rem.Mul(&rem, bigTwo) + half := rem.Cmp(&divisor) + if c.Rounding.ShouldAddOne(&d.Coeff, d.Negative, half) { + d.Coeff.Add(&d.Coeff, bigOne) + // The coefficient changed, so recompute num digits in + // setExponent. + nd = unknownNumDigits } } } - // 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, unknownNumDigits, res, int64(x.Exponent), int64(-y.Exponent), -adjust, diff) - quo.Negative = neg - d.Set(&quo) + res |= d.setExponent(c, nd, res, shift, -adjCoeffs, -adjExp10) + d.Reduce(d) // remove trailing zeros return c.goError(res) } @@ -567,7 +552,7 @@ func (c *Context) Sqrt(d, x *Decimal) (Condition, error) { d.Exponent += int32(e / 2) nc.Precision = c.Precision nc.Rounding = RoundHalfEven - d.Reduce(d) // Remove trailing zeros. + d.Reduce(d) // remove trailing zeros res := nc.round(d, d) return nc.goError(res) } diff --git a/decimal.go b/decimal.go index 03c8a88..49e9f0f 100644 --- a/decimal.go +++ b/decimal.go @@ -716,7 +716,7 @@ func (d *Decimal) Reduce(x *Decimal) (*Decimal, int) { d.Set(x) // Use a uint64 for the division if possible. - if d.Coeff.BitLen() <= 64 { + if d.Coeff.IsUint64() { i := d.Coeff.Uint64() for i >= 10000 && i%10000 == 0 { i /= 10000 diff --git a/gda_test.go b/gda_test.go index 177ff13..159dcf2 100644 --- a/gda_test.go +++ b/gda_test.go @@ -711,10 +711,6 @@ func (tc TestCase) PrintIgnore() { } var GDAignore = map[string]bool{ - // GDA says decNumber should skip these - "powx4302": true, - "powx4303": true, - // Invalid context "expx901": true, "expx902": true, @@ -775,21 +771,6 @@ var GDAignore = map[string]bool{ // TODO(mjibson): fix tests below - // log10(x) with large exponents, overflows - "logx0020": true, - "logx1146": true, - "logx1147": true, - "logx1156": true, - "logx1157": true, - "logx1176": true, - "logx1177": true, - - // log10(x) where x is 1.0 +/- some tiny epsilon. Our results are close but - // differ in the last places. - "logx1304": true, - "logx1309": true, - "logx1310": true, - // overflows to infinity "powx4125": true, "powx4145": true, @@ -819,90 +800,13 @@ var GDAignore = map[string]bool{ "addx61618": true, // extreme input range, but should work - "lnx0902": true, - "sqtx8137": true, - "sqtx8139": true, - "sqtx8145": true, - "sqtx8147": true, - "sqtx8149": true, - "sqtx8150": true, - "sqtx8151": true, - "sqtx8158": true, - "sqtx8165": true, - "sqtx8168": true, - "sqtx8174": true, - "sqtx8175": true, - "sqtx8179": true, - "sqtx8182": true, - "sqtx8185": true, - "sqtx8186": true, - "sqtx8195": true, - "sqtx8196": true, - "sqtx8197": true, - "sqtx8199": true, - "sqtx8204": true, - "sqtx8212": true, - "sqtx8213": true, - "sqtx8214": true, - "sqtx8218": true, - "sqtx8219": true, - "sqtx8323": true, - "sqtx8324": true, - "sqtx8331": true, - "sqtx8626": true, - "sqtx8627": true, - "sqtx8628": true, - "sqtx8631": true, - "sqtx8632": true, - "sqtx8633": true, - "sqtx8634": true, "sqtx8636": true, - "sqtx8637": true, - "sqtx8639": true, - "sqtx8640": true, - "sqtx8641": true, "sqtx8644": true, - "sqtx8645": true, "sqtx8646": true, "sqtx8647": true, "sqtx8648": true, - "sqtx8649": true, "sqtx8650": true, "sqtx8651": true, - "sqtx8652": true, - "sqtx8654": true, - - // tricky cases of underflow subnormals - "sqtx8700": true, - "sqtx8701": true, - "sqtx8702": true, - "sqtx8703": true, - "sqtx8704": true, - "sqtx8705": true, - "sqtx8706": true, - "sqtx8707": true, - "sqtx8708": true, - "sqtx8709": true, - "sqtx8710": true, - "sqtx8711": true, - "sqtx8712": true, - "sqtx8713": true, - "sqtx8714": true, - "sqtx8715": true, - "sqtx8716": true, - "sqtx8717": true, - "sqtx8718": true, - "sqtx8719": true, - "sqtx8720": true, - "sqtx8721": true, - "sqtx8722": true, - "sqtx8723": true, - "sqtx8724": true, - "sqtx8725": true, - "sqtx8726": true, - "sqtx8727": true, - "sqtx8728": true, - "sqtx8729": true, } var GDAignoreFlags = map[string]bool{ @@ -915,9 +819,4 @@ var GDAignoreFlags = map[string]bool{ "sqtx9039": true, "sqtx9040": true, "sqtx9045": true, - - // missing underflow, subnormal - "expx048": true, - "expx756": true, - "expx757": true, }