Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

apd: optimize Quo, remove per-digit long division #114

Merged
merged 5 commits into from
Jan 31, 2022
Merged

apd: optimize Quo, remove per-digit long division #114

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
da6c77f
apd: use IsUint64 instead of BitLen in Decimal.Reduce
nvanbenschoten Jan 31, 2022
909da8f
apd: optimize Quo, remove per-digit long-division
nvanbenschoten Jan 31, 2022
ae1b836
apd: remove looping when adjusting coefficients in Quo
nvanbenschoten Jan 31, 2022
c46294f
apd: unskip some GDA tests
nvanbenschoten Jan 31, 2022
d607123
apd: remove precision limit in Quo
nvanbenschoten Jan 31, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
151 changes: 68 additions & 83 deletions context.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(&dividend, 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(&dividend) <= 0 {
dividend.Sub(&dividend, &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(&dividend, 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(&dividend)
ndDivisor := NumDigits(&divisor)
ndDiff := ndDividend - ndDivisor
var tmpE BigInt
if ndDiff < 0 {
// numDigits(dividend) < numDigits(divisor), multiply dividend by 10^diff.
dividend.Mul(&dividend, 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(&dividend, 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(&dividend, tableExp10(adjExp10, &tmpE))

// Perform the division.
var rem BigInt
d.Coeff.QuoRem(&dividend, &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(&dividend, 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)
}

Expand Down Expand Up @@ -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)
}
Expand Down
2 changes: 1 addition & 1 deletion decimal.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
101 changes: 0 additions & 101 deletions gda_test.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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{
Expand All @@ -915,9 +819,4 @@ var GDAignoreFlags = map[string]bool{
"sqtx9039": true,
"sqtx9040": true,
"sqtx9045": true,

// missing underflow, subnormal
"expx048": true,
"expx756": true,
"expx757": true,
}