From 87eba72206a834c6822aed6aa9f92b0d21345e9a Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 27 May 2025 21:31:18 +0900 Subject: [PATCH 1/2] Use the same division/rounding logic in both BigDecimal_div and BigDecimal_div2 --- ext/bigdecimal/bigdecimal.c | 127 ++++++++++------------------- test/bigdecimal/test_bigdecimal.rb | 10 +-- 2 files changed, 46 insertions(+), 91 deletions(-) diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c index a707295..c18217e 100644 --- a/ext/bigdecimal/bigdecimal.c +++ b/ext/bigdecimal/bigdecimal.c @@ -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: @@ -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); @@ -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; @@ -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 = 2 * ((a_prec > b_prec) ? a_prec : b_prec); + 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; + GUARD_OBJ(res, NewZeroWrapNolimit(1, mx * VpBaseFig())); + VpDivd(cv, res, av, bv); + VpSetPrecLimit(pl); + VpLeftRound(cv, VpGetRoundMode(), ix); + return VpCheckGetValue(cv); } /* @@ -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; ind_r = 1; r->frac[0] = 0; diff --git a/test/bigdecimal/test_bigdecimal.rb b/test/bigdecimal/test_bigdecimal.rb index fab9622..cdebe42 100644 --- a/test/bigdecimal/test_bigdecimal.rb +++ b/test/bigdecimal/test_bigdecimal.rb @@ -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) 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) @@ -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.2722883438925926879095201027489267529117792091813217447000327237290151516073e-9") assert_equal(c, x / y, "[GH-220]") end @@ -1119,7 +1119,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 @@ -1129,7 +1129,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)) From b38d01521f062a79818dffedddb3eb7105bd8fea Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 27 May 2025 22:58:57 +0900 Subject: [PATCH 2/2] Reduce unnecessarily high division precision when precision is not specified Fixes exponential precision grouth. `max(divisor_prec, dividend_prec) + extra_prec` is enough for division precision. --- ext/bigdecimal/bigdecimal.c | 2 +- test/bigdecimal/test_bigdecimal.rb | 19 ++++++++++++++++++- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c index c18217e..954f627 100644 --- a/ext/bigdecimal/bigdecimal.c +++ b/ext/bigdecimal/bigdecimal.c @@ -2166,7 +2166,7 @@ BigDecimal_div2(VALUE self, VALUE b, VALUE n) 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 = 2 * ((a_prec > b_prec) ? a_prec : b_prec); + ix = ((a_prec > b_prec) ? a_prec : b_prec) + BIGDECIMAL_DOUBLE_FIGURES; if (2 * BIGDECIMAL_DOUBLE_FIGURES > ix) ix = 2 * BIGDECIMAL_DOUBLE_FIGURES; } diff --git a/test/bigdecimal/test_bigdecimal.rb b/test/bigdecimal/test_bigdecimal.rb index cdebe42..ac14a0e 100644 --- a/test/bigdecimal/test_bigdecimal.rb +++ b/test/bigdecimal/test_bigdecimal.rb @@ -982,7 +982,7 @@ def test_div def test_div_gh220 x = BigDecimal("1.0") y = BigDecimal("3672577333.6608990499165058135986328125") - c = BigDecimal("0.2722883438925926879095201027489267529117792091813217447000327237290151516073e-9") + c = BigDecimal("0.272288343892592687909520102748926752911779209181321745e-9") assert_equal(c, x / y, "[GH-220]") end @@ -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)