Skip to content

Commit 922fa0a

Browse files
authored
Fix VpDivd to fully use quotient array (#377)
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.
1 parent 4d7bb5b commit 922fa0a

File tree

5 files changed

+175
-19
lines changed

5 files changed

+175
-19
lines changed

.github/workflows/ci.yml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,7 @@ jobs:
4444
- { os: windows-latest , ruby: truffleruby-head }
4545
env:
4646
BIGDECIMAL_USE_DECDIG_UINT16_T: ${{ matrix.decdig_bits == 16 }}
47+
BIGDECIMAL_USE_VP_TEST_METHODS: true
4748

4849
steps:
4950
- uses: actions/checkout@v4

ext/bigdecimal/bigdecimal.c

Lines changed: 47 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -2158,11 +2158,10 @@ BigDecimal_div2(VALUE self, VALUE b, VALUE n)
21582158
ix = 2 * BIGDECIMAL_DOUBLE_FIGURES;
21592159
}
21602160

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

2164-
mx = bv.real->Prec + cv.real->MaxPrec - 1;
2165-
if (mx <= av.real->Prec) mx = av.real->Prec + 1;
2164+
mx = Max(av.real->Prec, bv.real->Prec + cv.real->MaxPrec - 1);
21662165
GUARD_OBJ(res, NewZeroWrapNolimit(1, mx * VpBaseFig()));
21672166
VpDivd(cv.real, res.real, av.real, bv.real);
21682167
VpSetPrecLimit(pl);
@@ -4215,6 +4214,35 @@ BigDecimal_literal(const char *str)
42154214

42164215
#define BIGDECIMAL_LITERAL(var, val) (BIGDECIMAL_ ## var = BigDecimal_literal(#val))
42174216

4217+
#ifdef BIGDECIMAL_USE_VP_TEST_METHODS
4218+
VALUE
4219+
BigDecimal_vpdivd(VALUE self, VALUE r, VALUE cprec) {
4220+
BDVALUE a,b,c,d;
4221+
size_t cn = NUM2INT(cprec);
4222+
a = GetBDValueMust(self);
4223+
b = GetBDValueMust(r);
4224+
size_t dn = Max(a.real->Prec, b.real->Prec + cn - 1);
4225+
c = NewZeroWrapLimited(1, cn * BASE_FIG);
4226+
d = NewZeroWrapLimited(1, dn * BASE_FIG);
4227+
VpDivd(c.real, d.real, a.real, b.real);
4228+
VpNmlz(c.real);
4229+
VpNmlz(d.real);
4230+
return rb_assoc_new(c.bigdecimal, d.bigdecimal);
4231+
}
4232+
4233+
VALUE
4234+
BigDecimal_vpmult(VALUE self, VALUE v) {
4235+
BDVALUE a,b,c;
4236+
a = GetBDValueMust(self);
4237+
b = GetBDValueMust(v);
4238+
size_t cn = a.real->Prec + b.real->Prec + 1;
4239+
c = NewZeroWrapLimited(1, cn * BASE_FIG);
4240+
VpMult(c.real, a.real, b.real);
4241+
VpNmlz(c.real);
4242+
return c.bigdecimal;
4243+
}
4244+
#endif /* BIGDECIMAL_USE_VP_TEST_METHODS */
4245+
42184246
/* Document-class: BigDecimal
42194247
* BigDecimal provides arbitrary-precision floating point decimal arithmetic.
42204248
*
@@ -4584,6 +4612,11 @@ Init_bigdecimal(void)
45844612
rb_define_method(rb_cBigDecimal, "truncate", BigDecimal_truncate, -1);
45854613
rb_define_method(rb_cBigDecimal, "_dump", BigDecimal_dump, -1);
45864614

4615+
#ifdef BIGDECIMAL_USE_VP_TEST_METHODS
4616+
rb_define_method(rb_cBigDecimal, "vpdivd", BigDecimal_vpdivd, 2);
4617+
rb_define_method(rb_cBigDecimal, "vpmult", BigDecimal_vpmult, 1);
4618+
#endif /* BIGDECIMAL_USE_VP_TEST_METHODS */
4619+
45874620
rb_mBigMath = rb_define_module("BigMath");
45884621
rb_define_singleton_method(rb_mBigMath, "exp", BigMath_s_exp, 2);
45894622
rb_define_singleton_method(rb_mBigMath, "log", BigMath_s_log, 2);
@@ -4645,7 +4678,6 @@ static int AddExponent(Real *a, SIGNED_VALUE n);
46454678
static DECDIG VpAddAbs(Real *a,Real *b,Real *c);
46464679
static DECDIG VpSubAbs(Real *a,Real *b,Real *c);
46474680
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);
4648-
static int VpNmlz(Real *a);
46494681
static void VpFormatSt(char *psz, size_t fFmt);
46504682
static int VpRdup(Real *m, size_t ind_m);
46514683

@@ -5977,6 +6009,10 @@ VpMult(Real *c, Real *a, Real *b)
59776009

59786010
/*
59796011
* c = a / b, remainder = r
6012+
* XXXX_YYYY_ZZZZ / 0001 = XXXX_YYYY_ZZZZ
6013+
* XXXX_YYYY_ZZZZ / 1111 = 000X_000Y_000Z
6014+
* 00XX_XXYY_YYZZ / 1000 = 0000_0XXX_XYYY
6015+
* 0001_0000_0000 / 9999 = 0000_0001_0001
59806016
*/
59816017
VP_EXPORT size_t
59826018
VpDivd(Real *c, Real *r, Real *a, Real *b)
@@ -6010,18 +6046,11 @@ VpDivd(Real *c, Real *r, Real *a, Real *b)
60106046
word_c = c->MaxPrec;
60116047
word_r = r->MaxPrec;
60126048

6013-
if (word_a >= word_r || word_b + word_c - 2 >= word_r) goto space_error;
6014-
6015-
ind_r = 1;
6016-
r->frac[0] = 0;
6017-
while (ind_r <= word_a) {
6018-
r->frac[ind_r] = a->frac[ind_r - 1];
6019-
++ind_r;
6020-
}
6021-
while (ind_r < word_r) r->frac[ind_r++] = 0;
6049+
if (word_a > word_r || word_b + word_c - 2 >= word_r) goto space_error;
60226050

6023-
ind_c = 0;
6024-
while (ind_c < word_c) c->frac[ind_c++] = 0;
6051+
for (i = 0; i < word_a; ++i) r->frac[i] = a->frac[i];
6052+
for (i = word_a; i < word_r; ++i) r->frac[i] = 0;
6053+
for (i = 0; i < word_c; ++i) c->frac[i] = 0;
60256054

60266055
/* initial procedure */
60276056
b1 = b1p1 = b->frac[0];
@@ -6037,7 +6066,7 @@ VpDivd(Real *c, Real *r, Real *a, Real *b)
60376066
/* */
60386067
/* loop start */
60396068
nLoop = Min(word_c, word_r);
6040-
ind_c = 1;
6069+
ind_c = 0;
60416070
while (ind_c < nLoop) {
60426071
if (r->frac[ind_c] == 0) {
60436072
++ind_c;
@@ -6141,13 +6170,12 @@ VpDivd(Real *c, Real *r, Real *a, Real *b)
61416170
c->Prec = word_c;
61426171
c->exponent = a->exponent;
61436172
VpSetSign(c, VpGetSign(a) * VpGetSign(b));
6144-
if (!AddExponent(c, 2)) return 0;
6173+
if (!AddExponent(c, 1)) return 0;
61456174
if (!AddExponent(c, -(b->exponent))) return 0;
61466175

61476176
VpNmlz(c); /* normalize c */
61486177
r->Prec = word_r;
61496178
r->exponent = a->exponent;
6150-
if (!AddExponent(r, 1)) return 0;
61516179
VpSetSign(r, VpGetSign(a));
61526180
VpNmlz(r); /* normalize r(remainder) */
61536181
goto Exit;

ext/bigdecimal/bigdecimal.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -220,6 +220,7 @@ VP_EXPORT size_t VpAsgn(Real *c, Real *a, int isw);
220220
VP_EXPORT size_t VpAddSub(Real *c,Real *a,Real *b,int operation);
221221
VP_EXPORT size_t VpMult(Real *c,Real *a,Real *b);
222222
VP_EXPORT size_t VpDivd(Real *c,Real *r,Real *a,Real *b);
223+
VP_EXPORT int VpNmlz(Real *a);
223224
VP_EXPORT int VpComp(Real *a,Real *b);
224225
VP_EXPORT ssize_t VpExponent10(Real *a);
225226
VP_EXPORT void VpSzMantissa(Real *a, char *buf, size_t bufsize);

ext/bigdecimal/extconf.rb

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,7 @@ def have_builtin_func(name, check_expr, opt = "", &b)
6060
end
6161

6262
$defs.push '-DBIGDECIMAL_USE_DECDIG_UINT16_T' if ENV['BIGDECIMAL_USE_DECDIG_UINT16_T'] == 'true'
63+
$defs.push '-DBIGDECIMAL_USE_VP_TEST_METHODS' if ENV['BIGDECIMAL_USE_VP_TEST_METHODS'] == 'true'
6364

6465
create_makefile('bigdecimal') {|mf|
6566
mf << "BIGDECIMAL_RB = #{bigdecimal_rb}\n"

test/bigdecimal/test_vp_operation.rb

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
# frozen_string_literal: false
2+
require_relative 'helper'
3+
require 'bigdecimal'
4+
5+
class TestVpOperation < Test::Unit::TestCase
6+
include TestBigDecimalBase
7+
8+
def setup
9+
super
10+
unless BigDecimal.instance_methods.include?(:vpdivd)
11+
# rake clean && BIGDECIMAL_USE_VP_TEST_METHODS=true rake compile
12+
omit 'Compile with BIGDECIMAL_USE_VP_TEST_METHODS=true to run this test'
13+
end
14+
end
15+
16+
def test_vpmult
17+
assert_equal(BigDecimal('121932631112635269'), BigDecimal('123456789').vpmult(BigDecimal('987654321')))
18+
assert_equal(BigDecimal('12193263.1112635269'), BigDecimal('123.456789').vpmult(BigDecimal('98765.4321')))
19+
x = 123**456
20+
y = 987**123
21+
assert_equal(BigDecimal("#{x * y}e-300"), BigDecimal("#{x}e-100").vpmult(BigDecimal("#{y}e-200")))
22+
end
23+
24+
def test_vpdivd
25+
# a[0] > b[0]
26+
# XXXX_YYYY_ZZZZ / 1111 #=> 000X_000Y_000Z
27+
x1 = BigDecimal('2' * BASE_FIG + '3' * BASE_FIG + '4' * BASE_FIG + '5' * BASE_FIG + '6' * BASE_FIG)
28+
y = BigDecimal('1' * BASE_FIG)
29+
d1 = BigDecimal("2e#{BASE_FIG * 4}")
30+
d2 = BigDecimal("3e#{BASE_FIG * 3}") + d1
31+
d3 = BigDecimal("4e#{BASE_FIG * 2}") + d2
32+
d4 = BigDecimal("5e#{BASE_FIG}") + d3
33+
d5 = BigDecimal(6) + d4
34+
assert_equal([d1, x1 - d1 * y], x1.vpdivd(y, 1))
35+
assert_equal([d2, x1 - d2 * y], x1.vpdivd(y, 2))
36+
assert_equal([d3, x1 - d3 * y], x1.vpdivd(y, 3))
37+
assert_equal([d4, x1 - d4 * y], x1.vpdivd(y, 4))
38+
assert_equal([d5, x1 - d5 * y], x1.vpdivd(y, 5))
39+
40+
# a[0] < b[0]
41+
# 00XX_XXYY_YYZZ_ZZ00 / 1111 #=> 0000_0X00_0Y00_0Z00
42+
shift = BASE_FIG / 2
43+
x2 = BigDecimal('2' * BASE_FIG + '3' * BASE_FIG + '4' * BASE_FIG + '5' * BASE_FIG + '6' * BASE_FIG + '0' * shift)
44+
d1 = BigDecimal("2e#{4 * BASE_FIG + shift}")
45+
d2 = BigDecimal("3e#{3 * BASE_FIG + shift}") + d1
46+
d3 = BigDecimal("4e#{2 * BASE_FIG + shift}") + d2
47+
d4 = BigDecimal("5e#{BASE_FIG + shift}") + d3
48+
d5 = BigDecimal("6e#{shift}") + d4
49+
assert_equal([0, x2], x2.vpdivd(y, 1))
50+
assert_equal([d1, x2 - d1 * y], x2.vpdivd(y, 2))
51+
assert_equal([d2, x2 - d2 * y], x2.vpdivd(y, 3))
52+
assert_equal([d3, x2 - d3 * y], x2.vpdivd(y, 4))
53+
assert_equal([d4, x2 - d4 * y], x2.vpdivd(y, 5))
54+
assert_equal([d5, x2 - d5 * y], x2.vpdivd(y, 6))
55+
end
56+
57+
def test_vpdivd_large_quotient_prec
58+
# 0001 / 0003 = 0000_3333_3333
59+
assert_equal([BigDecimal('0.' + '3' * BASE_FIG * 9), BigDecimal("1e-#{9 * BASE_FIG}")], BigDecimal(1).vpdivd(BigDecimal(3), 10))
60+
# 1000 / 0003 = 0333_3333_3333
61+
assert_equal([BigDecimal('3' * (BASE_FIG - 1) + '.' + '3' * BASE_FIG * 9), BigDecimal("1e-#{9 * BASE_FIG}")], BigDecimal(BASE / 10).vpdivd(BigDecimal(3), 10))
62+
end
63+
64+
def test_vpdivd_with_one
65+
x = BigDecimal('1234.2468000001234')
66+
assert_equal([BigDecimal('1234'), BigDecimal('0.2468000001234')], x.vpdivd(BigDecimal(1), 1))
67+
assert_equal([BigDecimal('+1234.2468'), BigDecimal('+0.1234e-9')], (+x).vpdivd(BigDecimal(+1), 2))
68+
assert_equal([BigDecimal('-1234.2468'), BigDecimal('+0.1234e-9')], (+x).vpdivd(BigDecimal(-1), 2))
69+
assert_equal([BigDecimal('-1234.2468'), BigDecimal('-0.1234e-9')], (-x).vpdivd(BigDecimal(+1), 2))
70+
assert_equal([BigDecimal('+1234.2468'), BigDecimal('-0.1234e-9')], (-x).vpdivd(BigDecimal(-1), 2))
71+
end
72+
73+
def test_vpdivd_precisions
74+
xs = [5, 10, 20, 40].map {|n| 123 ** n }
75+
ys = [5, 10, 20, 40].map {|n| 321 ** n }
76+
xs.product(ys).each do |x, y|
77+
[1, 2, 10, 20].each do |n|
78+
xn = (x.digits.size + BASE_FIG - 1) / BASE_FIG
79+
yn = (y.digits.size + BASE_FIG - 1) / BASE_FIG
80+
base = BASE ** (n - xn + yn - 1)
81+
div = BigDecimal((x * base / y).to_i) / base
82+
assert_equal([div, x - y * div], BigDecimal(x).vpdivd(y, n))
83+
end
84+
end
85+
end
86+
87+
def test_vpdivd_carry_borrow
88+
y_small = BASE / 7 * BASE ** 4
89+
y_large = (4 * BASE_FIG).times.map {|i| i % 9 + 1 }.join.to_i
90+
[y_large, y_small].each do |y|
91+
[0, 1, 2, BASE - 2, BASE - 1].repeated_permutation(4) do |a, b, c, d|
92+
x = y * (3 * BASE**4 + a * BASE**3 + b * BASE**2 + c * BASE + d) / BASE
93+
div = BigDecimal(x * BASE / y) / BASE
94+
mod = BigDecimal(x) - div * y
95+
assert_equal([div, mod], BigDecimal(x).vpdivd(BigDecimal(y), 5))
96+
end
97+
end
98+
end
99+
100+
def test_vpdivd_large_prec_divisor
101+
x = BigDecimal('2468.000000000000000000000000003')
102+
y1 = BigDecimal('1234.000000000000000000000000001')
103+
y2 = BigDecimal('1234.000000000000000000000000004')
104+
divy1_1 = BigDecimal(2)
105+
divy2_1 = BigDecimal(1)
106+
divy2_2 = BigDecimal('1.' + '9' * BASE_FIG)
107+
assert_equal([divy1_1, x - y1 * divy1_1], x.vpdivd(y1, 1))
108+
assert_equal([divy2_1, x - y2 * divy2_1], x.vpdivd(y2, 1))
109+
assert_equal([divy2_2, x - y2 * divy2_2], x.vpdivd(y2, 2))
110+
end
111+
112+
def test_vpdivd_intermediate_zero
113+
if BASE_FIG == 9
114+
x = BigDecimal('123456789.246913578000000000123456789')
115+
y = BigDecimal('123456789')
116+
assert_equal([BigDecimal('1.000000002000000000000000001'), BigDecimal(0)], x.vpdivd(y, 4))
117+
assert_equal([BigDecimal('1.000000000049999999'), BigDecimal('1e-18')], BigDecimal("2.000000000099999999").vpdivd(2, 3))
118+
else
119+
x = BigDecimal('1234.246800001234')
120+
y = BigDecimal('1234')
121+
assert_equal([BigDecimal('1.000200000001'), BigDecimal(0)], x.vpdivd(y, 4))
122+
assert_equal([BigDecimal('1.00000499'), BigDecimal('1e-8')], BigDecimal("2.00000999").vpdivd(2, 3))
123+
end
124+
end
125+
end

0 commit comments

Comments
 (0)