Skip to content

Fix precision loss in Newton-Raphson inverse #552

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

Open
tompng wants to merge 1 commit into
base: master
Choose a base branch
from
+84 −34
Open

Fix precision loss in Newton-Raphson inverse #552

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
74 changes: 54 additions & 20 deletions ext/bigdecimal/div.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,18 +1,45 @@
struct vp_settings {
size_t prec_limit;
unsigned short rounding_mode;
};

static VALUE
restore_vp_settings(VALUE saved_ptr)
{
struct vp_settings *saved = (struct vp_settings *)saved_ptr;
VpSetPrecLimit(saved->prec_limit);
VpSetRoundMode(saved->rounding_mode);
return Qnil;
}

// Calculate the inverse of x using the Newton-Raphson method.
// Assumes no precision limit and ROUND_HALF_UP. See VpDivdNewton.
static VALUE
newton_raphson_inverse(VALUE x, size_t prec) {
BDVALUE bdone = NewZeroWrap(1, 1);
VpSetOne(bdone.real);
VALUE one = bdone.bigdecimal;

// Initial approximation in 2 digits
// Initial approximation: 10^18 / (leading 9 digits of x), calculated in 64-bit integer.
// Truncating x to 9 digits and the quotient to an integer keep the relative error below 1.2e-8,
// so it has at least 7 correct digits.
const size_t initial_digits = 7;
const DECDIG_DBL base_sq = (DECDIG_DBL)BIGDECIMAL_BASE * BIGDECIMAL_BASE;
BDVALUE bdx = GetBDValueMust(x);
DECDIG_DBL x_lead = (DECDIG_DBL)bdx.real->frac[0] * BIGDECIMAL_BASE + (bdx.real->Prec >= 2 ? bdx.real->frac[1] : 0);
int shift = 0;
while (x_lead < base_sq / 10) {
x_lead *= 10;
shift++;
}
// Dividing base_sq - 1 instead of base_sq keeps the quotient below 10 * BASE,
// so that frac[0] of inv0 never reaches BASE.
DECDIG_DBL inv_lead = (base_sq - 1) / (x_lead / BIGDECIMAL_BASE);
for (int i = 0; i < shift; i++) inv_lead *= 10;
BDVALUE inv0 = NewZeroWrap(1, 2 * BIGDECIMAL_COMPONENT_FIGURES);
VpSetOne(inv0.real);
DECDIG_DBL numerator = (DECDIG_DBL)BIGDECIMAL_BASE * 100;
DECDIG_DBL denominator = (DECDIG_DBL)bdx.real->frac[0] * 100 + (DECDIG_DBL)(bdx.real->Prec >= 2 ? bdx.real->frac[1] : 0) * 100 / BIGDECIMAL_BASE;
inv0.real->frac[0] = (DECDIG)(numerator / denominator);
inv0.real->frac[1] = (DECDIG)((numerator % denominator) * (BIGDECIMAL_BASE / 100) / denominator * 100);
inv0.real->frac[0] = (DECDIG)(inv_lead / BIGDECIMAL_BASE);
inv0.real->frac[1] = (DECDIG)(inv_lead % BIGDECIMAL_BASE);
inv0.real->Prec = 2;
inv0.real->exponent = 1 - bdx.real->exponent;
VpNmlz(inv0.real);
Expand All @@ -22,14 +49,21 @@ newton_raphson_inverse(VALUE x, size_t prec) {
int bl = 1;
while (((size_t)1 << bl) < prec) bl++;

// Each iteration doubles the number of correct digits, and the rounding error of
// the previous iteration is squared into the new digits. Margin digits keep the
// squared error below the last digit, otherwise it can grow over iterations.
const size_t margin = 4;
for (int i = bl; i >= 0; i--) {
size_t n = (prec >> i) + 2;
size_t n = (prec >> i) + margin;
if (n > prec) n = prec;
// inv0 already has this precision. The last iteration is kept to round inv to prec digits.
if (n <= initial_digits && i > 0) continue;
// Newton-Raphson iteration: inv_next = inv + inv * (1 - x * inv)
// (1 - x * inv) is about 10^(-n/2), so calculating it in n/2 digits is enough for inv_next in n digits.
VALUE one_minus_x_inv = BigDecimal_sub2(
one,
BigDecimal_mult(BigDecimal_mult2(x, one, SIZET2NUM(n + 1)), inv),
SIZET2NUM(SIZET2NUM(n / 2))
SIZET2NUM(n / 2 + margin)
);
inv = BigDecimal_add2(
inv,
Expand Down Expand Up @@ -132,11 +166,15 @@ divmod_newton(VALUE x, VALUE y, VALUE *div_out, VALUE *mod_out) {
*mod_out = mod;
}

struct vp_divd_newton_args {
Real *c, *r, *a, *b;
};

static VALUE
VpDivdNewtonInner(VALUE args_ptr)
{
Real **args = (Real**)args_ptr;
Real *c = args[0], *r = args[1], *a = args[2], *b = args[3];
struct vp_divd_newton_args *args = (struct vp_divd_newton_args *)args_ptr;
Real *c = args->c, *r = args->r, *a = args->a, *b = args->b;
BDVALUE a2, b2, c2, r2;
VALUE div, mod, a2_frac = Qnil;
size_t div_prec = c->MaxPrec - 1;
Expand Down Expand Up @@ -174,19 +212,15 @@ VpDivdNewtonInner(VALUE args_ptr)
return Qnil;
}

static VALUE
ensure_restore_prec_limit(VALUE limit)
{
VpSetPrecLimit(NUM2SIZET(limit));
return Qnil;
}

// Newton-Raphson iteration converges from below, so a rounding mode that rounds toward zero
// (ROUND_DOWN, ROUND_FLOOR) adds error in the same direction at every iteration.
// Calculate with ROUND_HALF_UP and without precision limit, and restore them even if an exception is raised.
static void
VpDivdNewton(Real *c, Real *r, Real *a, Real *b)
{
Real *args[4] = {c, r, a, b};
size_t pl = VpGetPrecLimit();
struct vp_divd_newton_args args = {c, r, a, b};
struct vp_settings saved = {VpGetPrecLimit(), VpGetRoundMode()};
VpSetPrecLimit(0);
// Ensure restoring prec limit because some methods used in VpDivdNewtonInner may raise an exception
rb_ensure(VpDivdNewtonInner, (VALUE)args, ensure_restore_prec_limit, SIZET2NUM(pl));
VpSetRoundMode(VP_ROUND_HALF_UP);
rb_ensure(VpDivdNewtonInner, (VALUE)&args, restore_vp_settings, (VALUE)&saved);
}
44 changes: 30 additions & 14 deletions test/bigdecimal/test_vp_operation.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,21 +44,17 @@ def test_nttmult
end

def test_newton_inverse
# newton_raphson_inverse assumes ROUND_HALF_UP and no precision limit (the default of this test)
xs = [BigDecimal(3), BigDecimal('123e50'), BigDecimal('13' * 44), BigDecimal('17' * 45), BigDecimal('19' * 46)]
%i[up half_up down].each do |rounding_mode|
BigDecimal.save_rounding_mode do
BigDecimal.mode(BigDecimal::ROUND_MODE, rounding_mode)
[*1..32, 50, 100, 200, 300].each do |prec|
xs.each do |x|
inv = x.newton_raphson_inverse(prec)
assert_in_delta(1, x * inv, BigDecimal("1e#{1 - prec}"))

high_precision_inv = inv * (2 - x * inv)
expected_inv = high_precision_inv.mult(1, prec)
last_digit = BigDecimal("1e#{expected_inv.exponent - prec}")
assert_include([expected_inv - last_digit, expected_inv, expected_inv + last_digit], inv)
end
end
[*1..32, 50, 100, 200, 300].each do |prec|
xs.each do |x|
inv = x.newton_raphson_inverse(prec)
assert_in_delta(1, x * inv, BigDecimal("1e#{1 - prec}"))

high_precision_inv = inv * (2 - x * inv)
expected_inv = high_precision_inv.mult(1, prec)
last_digit = BigDecimal("1e#{expected_inv.exponent - prec}")
assert_include([expected_inv - last_digit, expected_inv, expected_inv + last_digit], inv)
end
end
end
Expand Down Expand Up @@ -88,6 +84,26 @@ def test_not_affected_by_limit
end
end

def test_not_affected_by_rounding_mode
x_int = 123**135
y_int = 135**123
mod_int = 111**111
x = BigDecimal(x_int)
y = BigDecimal(y_int)
mod = BigDecimal(mod_int)
z = BigDecimal(x_int * y_int + mod_int)
prec = (z.exponent - 1) / BASE_FIG - (y.exponent - 1) / BASE_FIG + 1
%i[up down ceiling floor half_even].each do |rounding_mode|
BigDecimal.save_rounding_mode do
mode = BigDecimal.mode(BigDecimal::ROUND_MODE, rounding_mode)
assert_equal([x, mod], z.vpdivd(y, prec))
assert_equal(mode, BigDecimal.mode(BigDecimal::ROUND_MODE))
assert_equal([x, mod], z.vpdivd_newton(y, prec))
assert_equal(mode, BigDecimal.mode(BigDecimal::ROUND_MODE))
end
end
end

def assert_vpdivd_equal(expected_divmod, x_y_n)
x, *args = x_y_n
assert_equal(expected_divmod, x.vpdivd(*args))
Expand Down