Skip to content

Fix VpDivd to fully use quotient array #377

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
tompng merged 1 commit into from
Jul 16, 2025
Merged

Fix VpDivd to fully use quotient array #377

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
1 change: 1 addition & 0 deletions .github/workflows/ci.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
66 changes: 47 additions & 19 deletions ext/bigdecimal/bigdecimal.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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.
*
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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];
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
1 change: 1 addition & 0 deletions ext/bigdecimal/bigdecimal.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
1 change: 1 addition & 0 deletions ext/bigdecimal/extconf.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
125 changes: 125 additions & 0 deletions test/bigdecimal/test_vp_operation.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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