From b9949823a701c7409593329e80e00a4fe98ff0fa Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 15 Jul 2025 20:01:07 +0900 Subject: [PATCH] Fix VpDivd to fully use quotient array VpDivd was using 1-index in quotient and remainder digits array. c->frac[0] and r->frac[0] was always zero. Fix to use 0-index. Required MaxPrec of remainder reduces. Add Test for VpDivd expected calculation. --- .github/workflows/ci.yml | 1 + ext/bigdecimal/bigdecimal.c | 66 ++++++++++---- ext/bigdecimal/bigdecimal.h | 1 + ext/bigdecimal/extconf.rb | 1 + test/bigdecimal/test_vp_operation.rb | 125 +++++++++++++++++++++++++++ 5 files changed, 175 insertions(+), 19 deletions(-) create mode 100644 test/bigdecimal/test_vp_operation.rb diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index da88a17b..1b935355 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -44,6 +44,7 @@ jobs: - { os: windows-latest , ruby: truffleruby-head } env: BIGDECIMAL_USE_DECDIG_UINT16_T: ${{ matrix.decdig_bits == 16 }} + BIGDECIMAL_USE_VP_TEST_METHODS: true steps: - uses: actions/checkout@v4 diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c index 93c82068..bce7f830 100644 --- a/ext/bigdecimal/bigdecimal.c +++ b/ext/bigdecimal/bigdecimal.c @@ -2158,11 +2158,10 @@ BigDecimal_div2(VALUE self, VALUE b, VALUE n) ix = 2 * BIGDECIMAL_DOUBLE_FIGURES; } - // VpDivd needs 2 extra DECDIGs. One more is needed for rounding. - GUARD_OBJ(cv, NewZeroWrapLimited(1, ix + 3 * VpBaseFig())); + // VpDivd needs 2 extra DECDIGs for rounding. + GUARD_OBJ(cv, NewZeroWrapLimited(1, ix + 2 * VpBaseFig())); - mx = bv.real->Prec + cv.real->MaxPrec - 1; - if (mx <= av.real->Prec) mx = av.real->Prec + 1; + mx = Max(av.real->Prec, bv.real->Prec + cv.real->MaxPrec - 1); GUARD_OBJ(res, NewZeroWrapNolimit(1, mx * VpBaseFig())); VpDivd(cv.real, res.real, av.real, bv.real); VpSetPrecLimit(pl); @@ -4215,6 +4214,35 @@ BigDecimal_literal(const char *str) #define BIGDECIMAL_LITERAL(var, val) (BIGDECIMAL_ ## var = BigDecimal_literal(#val)) +#ifdef BIGDECIMAL_USE_VP_TEST_METHODS +VALUE +BigDecimal_vpdivd(VALUE self, VALUE r, VALUE cprec) { + BDVALUE a,b,c,d; + size_t cn = NUM2INT(cprec); + a = GetBDValueMust(self); + b = GetBDValueMust(r); + size_t dn = Max(a.real->Prec, b.real->Prec + cn - 1); + c = NewZeroWrapLimited(1, cn * BASE_FIG); + d = NewZeroWrapLimited(1, dn * BASE_FIG); + VpDivd(c.real, d.real, a.real, b.real); + VpNmlz(c.real); + VpNmlz(d.real); + return rb_assoc_new(c.bigdecimal, d.bigdecimal); +} + +VALUE +BigDecimal_vpmult(VALUE self, VALUE v) { + BDVALUE a,b,c; + a = GetBDValueMust(self); + b = GetBDValueMust(v); + size_t cn = a.real->Prec + b.real->Prec + 1; + c = NewZeroWrapLimited(1, cn * BASE_FIG); + VpMult(c.real, a.real, b.real); + VpNmlz(c.real); + return c.bigdecimal; +} +#endif /* BIGDECIMAL_USE_VP_TEST_METHODS */ + /* Document-class: BigDecimal * BigDecimal provides arbitrary-precision floating point decimal arithmetic. * @@ -4584,6 +4612,11 @@ Init_bigdecimal(void) rb_define_method(rb_cBigDecimal, "truncate", BigDecimal_truncate, -1); rb_define_method(rb_cBigDecimal, "_dump", BigDecimal_dump, -1); +#ifdef BIGDECIMAL_USE_VP_TEST_METHODS + rb_define_method(rb_cBigDecimal, "vpdivd", BigDecimal_vpdivd, 2); + rb_define_method(rb_cBigDecimal, "vpmult", BigDecimal_vpmult, 1); +#endif /* BIGDECIMAL_USE_VP_TEST_METHODS */ + rb_mBigMath = rb_define_module("BigMath"); rb_define_singleton_method(rb_mBigMath, "exp", BigMath_s_exp, 2); rb_define_singleton_method(rb_mBigMath, "log", BigMath_s_log, 2); @@ -4645,7 +4678,6 @@ static int AddExponent(Real *a, SIGNED_VALUE n); static DECDIG VpAddAbs(Real *a,Real *b,Real *c); static DECDIG VpSubAbs(Real *a,Real *b,Real *c); static size_t VpSetPTR(Real *a, Real *b, Real *c, size_t *a_pos, size_t *b_pos, size_t *c_pos, DECDIG *av, DECDIG *bv); -static int VpNmlz(Real *a); static void VpFormatSt(char *psz, size_t fFmt); static int VpRdup(Real *m, size_t ind_m); @@ -5977,6 +6009,10 @@ VpMult(Real *c, Real *a, Real *b) /* * c = a / b, remainder = r + * XXXX_YYYY_ZZZZ / 0001 = XXXX_YYYY_ZZZZ + * XXXX_YYYY_ZZZZ / 1111 = 000X_000Y_000Z + * 00XX_XXYY_YYZZ / 1000 = 0000_0XXX_XYYY + * 0001_0000_0000 / 9999 = 0000_0001_0001 */ VP_EXPORT size_t VpDivd(Real *c, Real *r, Real *a, Real *b) @@ -6010,18 +6046,11 @@ VpDivd(Real *c, Real *r, Real *a, Real *b) word_c = c->MaxPrec; word_r = r->MaxPrec; - if (word_a >= word_r || word_b + word_c - 2 >= word_r) goto space_error; - - ind_r = 1; - r->frac[0] = 0; - while (ind_r <= word_a) { - r->frac[ind_r] = a->frac[ind_r - 1]; - ++ind_r; - } - while (ind_r < word_r) r->frac[ind_r++] = 0; + if (word_a > word_r || word_b + word_c - 2 >= word_r) goto space_error; - ind_c = 0; - while (ind_c < word_c) c->frac[ind_c++] = 0; + for (i = 0; i < word_a; ++i) r->frac[i] = a->frac[i]; + for (i = word_a; i < word_r; ++i) r->frac[i] = 0; + for (i = 0; i < word_c; ++i) c->frac[i] = 0; /* initial procedure */ b1 = b1p1 = b->frac[0]; @@ -6037,7 +6066,7 @@ VpDivd(Real *c, Real *r, Real *a, Real *b) /* */ /* loop start */ nLoop = Min(word_c, word_r); - ind_c = 1; + ind_c = 0; while (ind_c < nLoop) { if (r->frac[ind_c] == 0) { ++ind_c; @@ -6141,13 +6170,12 @@ VpDivd(Real *c, Real *r, Real *a, Real *b) c->Prec = word_c; c->exponent = a->exponent; VpSetSign(c, VpGetSign(a) * VpGetSign(b)); - if (!AddExponent(c, 2)) return 0; + if (!AddExponent(c, 1)) return 0; if (!AddExponent(c, -(b->exponent))) return 0; VpNmlz(c); /* normalize c */ r->Prec = word_r; r->exponent = a->exponent; - if (!AddExponent(r, 1)) return 0; VpSetSign(r, VpGetSign(a)); VpNmlz(r); /* normalize r(remainder) */ goto Exit; diff --git a/ext/bigdecimal/bigdecimal.h b/ext/bigdecimal/bigdecimal.h index 8e21c1c9..8a007823 100644 --- a/ext/bigdecimal/bigdecimal.h +++ b/ext/bigdecimal/bigdecimal.h @@ -220,6 +220,7 @@ VP_EXPORT size_t VpAsgn(Real *c, Real *a, int isw); VP_EXPORT size_t VpAddSub(Real *c,Real *a,Real *b,int operation); VP_EXPORT size_t VpMult(Real *c,Real *a,Real *b); VP_EXPORT size_t VpDivd(Real *c,Real *r,Real *a,Real *b); +VP_EXPORT int VpNmlz(Real *a); VP_EXPORT int VpComp(Real *a,Real *b); VP_EXPORT ssize_t VpExponent10(Real *a); VP_EXPORT void VpSzMantissa(Real *a, char *buf, size_t bufsize); diff --git a/ext/bigdecimal/extconf.rb b/ext/bigdecimal/extconf.rb index 72f9843f..a4825103 100644 --- a/ext/bigdecimal/extconf.rb +++ b/ext/bigdecimal/extconf.rb @@ -60,6 +60,7 @@ def have_builtin_func(name, check_expr, opt = "", &b) end $defs.push '-DBIGDECIMAL_USE_DECDIG_UINT16_T' if ENV['BIGDECIMAL_USE_DECDIG_UINT16_T'] == 'true' +$defs.push '-DBIGDECIMAL_USE_VP_TEST_METHODS' if ENV['BIGDECIMAL_USE_VP_TEST_METHODS'] == 'true' create_makefile('bigdecimal') {|mf| mf << "BIGDECIMAL_RB = #{bigdecimal_rb}\n" diff --git a/test/bigdecimal/test_vp_operation.rb b/test/bigdecimal/test_vp_operation.rb new file mode 100644 index 00000000..1ce3a034 --- /dev/null +++ b/test/bigdecimal/test_vp_operation.rb @@ -0,0 +1,125 @@ +# frozen_string_literal: false +require_relative 'helper' +require 'bigdecimal' + +class TestVpOperation < Test::Unit::TestCase + include TestBigDecimalBase + + def setup + super + unless BigDecimal.instance_methods.include?(:vpdivd) + # rake clean && BIGDECIMAL_USE_VP_TEST_METHODS=true rake compile + omit 'Compile with BIGDECIMAL_USE_VP_TEST_METHODS=true to run this test' + end + end + + def test_vpmult + assert_equal(BigDecimal('121932631112635269'), BigDecimal('123456789').vpmult(BigDecimal('987654321'))) + assert_equal(BigDecimal('12193263.1112635269'), BigDecimal('123.456789').vpmult(BigDecimal('98765.4321'))) + x = 123**456 + y = 987**123 + assert_equal(BigDecimal("#{x * y}e-300"), BigDecimal("#{x}e-100").vpmult(BigDecimal("#{y}e-200"))) + end + + def test_vpdivd + # a[0] > b[0] + # XXXX_YYYY_ZZZZ / 1111 #=> 000X_000Y_000Z + x1 = BigDecimal('2' * BASE_FIG + '3' * BASE_FIG + '4' * BASE_FIG + '5' * BASE_FIG + '6' * BASE_FIG) + y = BigDecimal('1' * BASE_FIG) + d1 = BigDecimal("2e#{BASE_FIG * 4}") + d2 = BigDecimal("3e#{BASE_FIG * 3}") + d1 + d3 = BigDecimal("4e#{BASE_FIG * 2}") + d2 + d4 = BigDecimal("5e#{BASE_FIG}") + d3 + d5 = BigDecimal(6) + d4 + assert_equal([d1, x1 - d1 * y], x1.vpdivd(y, 1)) + assert_equal([d2, x1 - d2 * y], x1.vpdivd(y, 2)) + assert_equal([d3, x1 - d3 * y], x1.vpdivd(y, 3)) + assert_equal([d4, x1 - d4 * y], x1.vpdivd(y, 4)) + assert_equal([d5, x1 - d5 * y], x1.vpdivd(y, 5)) + + # a[0] < b[0] + # 00XX_XXYY_YYZZ_ZZ00 / 1111 #=> 0000_0X00_0Y00_0Z00 + shift = BASE_FIG / 2 + x2 = BigDecimal('2' * BASE_FIG + '3' * BASE_FIG + '4' * BASE_FIG + '5' * BASE_FIG + '6' * BASE_FIG + '0' * shift) + d1 = BigDecimal("2e#{4 * BASE_FIG + shift}") + d2 = BigDecimal("3e#{3 * BASE_FIG + shift}") + d1 + d3 = BigDecimal("4e#{2 * BASE_FIG + shift}") + d2 + d4 = BigDecimal("5e#{BASE_FIG + shift}") + d3 + d5 = BigDecimal("6e#{shift}") + d4 + assert_equal([0, x2], x2.vpdivd(y, 1)) + assert_equal([d1, x2 - d1 * y], x2.vpdivd(y, 2)) + assert_equal([d2, x2 - d2 * y], x2.vpdivd(y, 3)) + assert_equal([d3, x2 - d3 * y], x2.vpdivd(y, 4)) + assert_equal([d4, x2 - d4 * y], x2.vpdivd(y, 5)) + assert_equal([d5, x2 - d5 * y], x2.vpdivd(y, 6)) + end + + def test_vpdivd_large_quotient_prec + # 0001 / 0003 = 0000_3333_3333 + assert_equal([BigDecimal('0.' + '3' * BASE_FIG * 9), BigDecimal("1e-#{9 * BASE_FIG}")], BigDecimal(1).vpdivd(BigDecimal(3), 10)) + # 1000 / 0003 = 0333_3333_3333 + assert_equal([BigDecimal('3' * (BASE_FIG - 1) + '.' + '3' * BASE_FIG * 9), BigDecimal("1e-#{9 * BASE_FIG}")], BigDecimal(BASE / 10).vpdivd(BigDecimal(3), 10)) + end + + def test_vpdivd_with_one + x = BigDecimal('1234.2468000001234') + assert_equal([BigDecimal('1234'), BigDecimal('0.2468000001234')], x.vpdivd(BigDecimal(1), 1)) + assert_equal([BigDecimal('+1234.2468'), BigDecimal('+0.1234e-9')], (+x).vpdivd(BigDecimal(+1), 2)) + assert_equal([BigDecimal('-1234.2468'), BigDecimal('+0.1234e-9')], (+x).vpdivd(BigDecimal(-1), 2)) + assert_equal([BigDecimal('-1234.2468'), BigDecimal('-0.1234e-9')], (-x).vpdivd(BigDecimal(+1), 2)) + assert_equal([BigDecimal('+1234.2468'), BigDecimal('-0.1234e-9')], (-x).vpdivd(BigDecimal(-1), 2)) + end + + def test_vpdivd_precisions + xs = [5, 10, 20, 40].map {|n| 123 ** n } + ys = [5, 10, 20, 40].map {|n| 321 ** n } + xs.product(ys).each do |x, y| + [1, 2, 10, 20].each do |n| + xn = (x.digits.size + BASE_FIG - 1) / BASE_FIG + yn = (y.digits.size + BASE_FIG - 1) / BASE_FIG + base = BASE ** (n - xn + yn - 1) + div = BigDecimal((x * base / y).to_i) / base + assert_equal([div, x - y * div], BigDecimal(x).vpdivd(y, n)) + end + end + end + + def test_vpdivd_carry_borrow + y_small = BASE / 7 * BASE ** 4 + y_large = (4 * BASE_FIG).times.map {|i| i % 9 + 1 }.join.to_i + [y_large, y_small].each do |y| + [0, 1, 2, BASE - 2, BASE - 1].repeated_permutation(4) do |a, b, c, d| + x = y * (3 * BASE**4 + a * BASE**3 + b * BASE**2 + c * BASE + d) / BASE + div = BigDecimal(x * BASE / y) / BASE + mod = BigDecimal(x) - div * y + assert_equal([div, mod], BigDecimal(x).vpdivd(BigDecimal(y), 5)) + end + end + end + + def test_vpdivd_large_prec_divisor + x = BigDecimal('2468.000000000000000000000000003') + y1 = BigDecimal('1234.000000000000000000000000001') + y2 = BigDecimal('1234.000000000000000000000000004') + divy1_1 = BigDecimal(2) + divy2_1 = BigDecimal(1) + divy2_2 = BigDecimal('1.' + '9' * BASE_FIG) + assert_equal([divy1_1, x - y1 * divy1_1], x.vpdivd(y1, 1)) + assert_equal([divy2_1, x - y2 * divy2_1], x.vpdivd(y2, 1)) + assert_equal([divy2_2, x - y2 * divy2_2], x.vpdivd(y2, 2)) + end + + def test_vpdivd_intermediate_zero + if BASE_FIG == 9 + x = BigDecimal('123456789.246913578000000000123456789') + y = BigDecimal('123456789') + assert_equal([BigDecimal('1.000000002000000000000000001'), BigDecimal(0)], x.vpdivd(y, 4)) + assert_equal([BigDecimal('1.000000000049999999'), BigDecimal('1e-18')], BigDecimal("2.000000000099999999").vpdivd(2, 3)) + else + x = BigDecimal('1234.246800001234') + y = BigDecimal('1234') + assert_equal([BigDecimal('1.000200000001'), BigDecimal(0)], x.vpdivd(y, 4)) + assert_equal([BigDecimal('1.00000499'), BigDecimal('1e-8')], BigDecimal("2.00000999").vpdivd(2, 3)) + end + end +end