Skip to content

Integrate BigDecimal_div and BigDecimal_div2 #329

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

Merged
merged 2 commits into from
May 29, 2025
Merged

Integrate BigDecimal_div and BigDecimal_div2 #329

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
87eba72
Use the same division/rounding logic in both BigDecimal_div and BigDe…
tompng May 27, 2025
b38d015
Reduce unnecessarily high division precision when precision is not sp…
tompng May 27, 2025
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
127 changes: 41 additions & 86 deletions ext/bigdecimal/bigdecimal.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1818,55 +1818,6 @@ BigDecimal_mult(VALUE self, VALUE r)
return VpCheckGetValue(c);
}

static VALUE
BigDecimal_divide(VALUE self, VALUE r, Real **c, Real **res, Real **div)
/* For c = self.div(r): with round operation */
{
ENTER(5);
Real *a, *b;
ssize_t a_prec, b_prec;
size_t mx;

TypedData_Get_Struct(self, Real, &BigDecimal_data_type, a);
SAVE(a);

VALUE rr = r;
if (is_kind_of_BigDecimal(rr)) {
/* do nothing */
}
else if (RB_INTEGER_TYPE_P(r)) {
rr = rb_inum_convert_to_BigDecimal(r, 0, true);
}
else if (RB_TYPE_P(r, T_FLOAT)) {
rr = rb_float_convert_to_BigDecimal(r, 0, true);
}
else if (RB_TYPE_P(r, T_RATIONAL)) {
rr = rb_rational_convert_to_BigDecimal(r, a->Prec*BASE_FIG, true);
}

if (!is_kind_of_BigDecimal(rr)) {
return DoSomeOne(self, r, '/');
}

TypedData_Get_Struct(rr, Real, &BigDecimal_data_type, b);
SAVE(b);
*div = b;

BigDecimal_count_precision_and_scale(self, &a_prec, NULL);
BigDecimal_count_precision_and_scale(rr, &b_prec, NULL);
mx = (a_prec > b_prec) ? a_prec : b_prec;
mx *= 2;

if (2*BIGDECIMAL_DOUBLE_FIGURES > mx)
mx = 2*BIGDECIMAL_DOUBLE_FIGURES;

GUARD_OBJ((*c), NewZeroWrapNolimit(1, mx + 2*BASE_FIG));
GUARD_OBJ((*res), NewZeroWrapNolimit(1, (mx + 1)*2 + 2*BASE_FIG));
VpDivd(*c, *res, a, b);

return Qnil;
}

static VALUE BigDecimal_DoDivmod(VALUE self, VALUE r, Real **div, Real **mod);

/* call-seq:
Expand All @@ -1884,20 +1835,15 @@ static VALUE
BigDecimal_div(VALUE self, VALUE r)
/* For c = self/r: with round operation */
{
ENTER(5);
Real *c=NULL, *res=NULL, *div = NULL;
r = BigDecimal_divide(self, r, &c, &res, &div);
if (!NIL_P(r)) return r; /* coerced by other */
SAVE(c); SAVE(res); SAVE(div);
/* a/b = c + r/b */
/* c xxxxx
r 00000yyyyy ==> (y/b)*BASE >= HALF_BASE
*/
/* Round */
if (VpHasVal(div)) { /* frac[0] must be zero for NaN,INF,Zero */
VpInternalRound(c, 0, c->frac[c->Prec-1], (DECDIG)(VpBaseVal() * (DECDIG_DBL)res->frac[0] / div->frac[0]));
if (
!is_kind_of_BigDecimal(r) &&
!RB_INTEGER_TYPE_P(r) &&
!RB_TYPE_P(r, T_FLOAT) &&
!RB_TYPE_P(r, T_RATIONAL)
) {
return DoSomeOne(self, r, '/');
}
return VpCheckGetValue(c);
return BigDecimal_div2(self, r, INT2FIX(0));
}

static VALUE BigDecimal_round(int argc, VALUE *argv, VALUE self);
Expand Down Expand Up @@ -2187,6 +2133,9 @@ BigDecimal_div2(VALUE self, VALUE b, VALUE n)
{
ENTER(5);
SIGNED_VALUE ix;
Real *res = NULL;
Real *av = NULL, *bv = NULL, *cv = NULL;
size_t mx, pl;

if (NIL_P(n)) { /* div in Float sense */
Real *div = NULL;
Expand All @@ -2199,33 +2148,39 @@ BigDecimal_div2(VALUE self, VALUE b, VALUE n)

/* div in BigDecimal sense */
ix = check_int_precision(n);
if (ix == 0) {
return BigDecimal_div(self, b);
}
else {
Real *res = NULL;
Real *av = NULL, *bv = NULL, *cv = NULL;
size_t mx = ix + VpBaseFig()*2;
size_t b_prec = ix;
size_t pl = VpSetPrecLimit(0);

GUARD_OBJ(cv, NewZeroWrapLimited(1, mx + VpBaseFig()));
GUARD_OBJ(av, GetVpValue(self, 1));

pl = VpSetPrecLimit(0);

GUARD_OBJ(av, GetVpValue(self, 1));
if (RB_FLOAT_TYPE_P(b) && ix > BIGDECIMAL_DOUBLE_FIGURES) {
/* TODO: I want to refactor this precision control for a float value later
* by introducing an implicit conversion function instead of
* GetVpValueWithPrec. */
if (RB_FLOAT_TYPE_P(b) && b_prec > BIGDECIMAL_DOUBLE_FIGURES) {
b_prec = BIGDECIMAL_DOUBLE_FIGURES;
}
GUARD_OBJ(bv, GetVpValueWithPrec(b, b_prec, 1));
mx = av->Prec + bv->Prec + 2;
if (mx <= cv->MaxPrec) mx = cv->MaxPrec + 1;
GUARD_OBJ(res, NewZeroWrapNolimit(1, (mx * 2 + 2)*VpBaseFig()));
VpDivd(cv, res, av, bv);
VpSetPrecLimit(pl);
VpLeftRound(cv, VpGetRoundMode(), ix);
return VpCheckGetValue(cv);
GUARD_OBJ(bv, GetVpValueWithPrec(b, BIGDECIMAL_DOUBLE_FIGURES, 1));
}
else {
GUARD_OBJ(bv, GetVpValueWithPrec(b, ix, 1));
}

if (ix == 0) {
ssize_t a_prec, b_prec;
BigDecimal_count_precision_and_scale(av->obj, &a_prec, NULL);
BigDecimal_count_precision_and_scale(bv->obj, &b_prec, NULL);
ix = ((a_prec > b_prec) ? a_prec : b_prec) + BIGDECIMAL_DOUBLE_FIGURES;
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

max(a_prec, b_prec) + extra_prec is enough. Fixes exponential precision growth.

if (2 * BIGDECIMAL_DOUBLE_FIGURES > ix)
ix = 2 * BIGDECIMAL_DOUBLE_FIGURES;
}

// VpDivd needs 2 extra DECDIGs. One more is needed for rounding.
GUARD_OBJ(cv, NewZeroWrapLimited(1, ix + 3 * VpBaseFig()));

mx = bv->Prec + cv->MaxPrec - 1;
if (mx <= av->Prec) mx = av->Prec + 1;
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Correct allocation size calculation of res.
VpDivd requires word_a < word_r && word_b + word_c - 2 < word_r

GUARD_OBJ(res, NewZeroWrapNolimit(1, mx * VpBaseFig()));
VpDivd(cv, res, av, bv);
VpSetPrecLimit(pl);
VpLeftRound(cv, VpGetRoundMode(), ix);
return VpCheckGetValue(cv);
}

/*
Expand Down Expand Up @@ -6164,7 +6119,7 @@ VpDivd(Real *c, Real *r, Real *a, Real *b)
word_c = c->MaxPrec;
word_r = r->MaxPrec;

if (word_a >= word_r) goto space_error;
if (word_a >= word_r || word_b + word_c - 2 >= word_r) goto space_error;
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are two more goto space_error in this function.
One is on L6218

ind_r = ind_c + ind_b;
if (ind_r >= word_r) goto space_error;

And another in L6253

ind_r = b->Prec + ind_c - 1; // L6242
ind_r = b->Prec + ind_c; // L6250 (ind_c + 1 < word_c is ensured in L6247)
if (ind_r >= word_r) goto space_error;

Both has the same requirement: (word_b-1) + (word_c-1) < word_r


ind_r = 1;
r->frac[0] = 0;
Expand Down
27 changes: 22 additions & 5 deletions test/bigdecimal/test_bigdecimal.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -510,10 +510,10 @@ def test_round_up

BigDecimal.mode(BigDecimal::ROUND_MODE, BigDecimal::ROUND_HALF_DOWN)
assert_operator(n4, :>, n4 / 9 * 9)
assert_operator(n5, :>, n5 / 9 * 9)
assert_operator(n5, :<, n5 / 9 * 9)
Copy link
Member Author

@tompng tompng May 27, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Half down rules are:

123.5 → 123 (half down)
123.50001 → 124 (round up because it's larger than half)
0.55550 → 0.555
0.55551 → 0.556
0.55555555.. → 0.556 (Because it's larger than half)

In master branch, rounding is inconsistent between a / b and a.div(b, prec).

BigDecimal.mode(BigDecimal::ROUND_MODE, BigDecimal::ROUND_HALF_DOWN)

BigDecimal(5)/9
#=> 0.555555555555555555555555555555555555e0 (ends with 5)
BigDecimal(5).div(9, 36)
#=> 0.555555555555555555555555555555555556e0 (ends with 6)

assert_operator(n6, :<, n6 / 9 * 9)
assert_operator(m4, :<, m4 / 9 * 9)
assert_operator(m5, :<, m5 / 9 * 9)
assert_operator(m5, :>, m5 / 9 * 9)
assert_operator(m6, :>, m6 / 9 * 9)
assert_equal(2, n2h.round)
assert_equal(3, n3h.round)
Expand Down Expand Up @@ -982,7 +982,7 @@ def test_div
def test_div_gh220
x = BigDecimal("1.0")
y = BigDecimal("3672577333.6608990499165058135986328125")
c = BigDecimal("0.272288343892592687909520102748926752911779209181321744700032723729015151607289998e-9")
c = BigDecimal("0.272288343892592687909520102748926752911779209181321745e-9")
assert_equal(c, x / y, "[GH-220]")
end

Expand All @@ -995,6 +995,23 @@ def test_div_precision
"(101/0.9163472602589686).precision >= (0.9163472602589686).precision #{bug13754}")
end

def test_div_various_precisions
a_precs = [5, 20, 70]
b_precs = [*5..80]
exponents = [-5, 0, 5]
a_precs.product(exponents, b_precs, exponents).each do |prec_a, ex_a, prec_b, ex_b|
a = BigDecimal('7.' + '1' * (prec_a - 1) + "e#{ex_a}")
b = BigDecimal('3.' + '1' * (prec_b - 1) + "e#{ex_b}")
c = a / b
max = [prec_a, prec_b, BigDecimal.double_fig].max
# Precision must be enough and not too large
precision_min = max + BigDecimal.double_fig / 2
precision_max = max + 2 * BigDecimal.double_fig
assert_includes(precision_min..precision_max, c.n_significant_digits)
assert_in_delta(a, c * b, a * 10**(1 - precision_min))
end
end

def test_div_with_float
assert_kind_of(BigDecimal, BigDecimal("3") / 1.5)
assert_equal(BigDecimal("0.5"), BigDecimal(1) / 2.0)
Expand Down Expand Up @@ -1119,7 +1136,7 @@ def test_div_bigdecimal_with_float_and_precision
def test_quo_without_prec
x = BigDecimal(5)
y = BigDecimal(229)
assert_equal(BigDecimal("0.021834061135371179039301310043668122"), x.quo(y))
assert_equal(BigDecimal("0.021834061135371179039301310043668"), x.quo(y))
end

def test_quo_with_prec
Expand All @@ -1129,7 +1146,7 @@ def test_quo_with_prec

x = BigDecimal(5)
y = BigDecimal(229)
assert_equal(BigDecimal("0.021834061135371179039301310043668122"), x.quo(y, 0))
assert_equal(BigDecimal("0.021834061135371179039301310043668"), x.quo(y, 0))
assert_equal(BigDecimal("0.022"), x.quo(y, 2))
assert_equal(BigDecimal("0.0218"), x.quo(y, 3))
assert_equal(BigDecimal("0.0218341"), x.quo(y, 6))
Expand Down