From ed3ea560c4b99b337b093c5919864250f5172adc Mon Sep 17 00:00:00 2001 From: tompng Date: Mon, 19 May 2025 21:04:51 +0900 Subject: [PATCH 1/2] Implement BigDecimal#_decimal_shift for internal use --- ext/bigdecimal/bigdecimal.c | 60 +++++++++++++++++++++++++++++- test/bigdecimal/test_bigdecimal.rb | 17 +++++++++ 2 files changed, 76 insertions(+), 1 deletion(-) diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c index 70cc4196..b63c6735 100644 --- a/ext/bigdecimal/bigdecimal.c +++ b/ext/bigdecimal/bigdecimal.c @@ -207,6 +207,7 @@ rbd_allocate_struct_zero(int sign, size_t const digits) static unsigned short VpGetException(void); static void VpSetException(unsigned short f); static void VpCheckException(Real *p, bool always); +static int AddExponent(Real *a, SIGNED_VALUE n); static VALUE CheckGetValue(BDVALUE v); static void VpInternalRound(Real *c, size_t ixDigit, DECDIG vPrev, DECDIG v); static int VpLimitRound(Real *c, size_t ixDigit); @@ -2453,6 +2454,63 @@ BigDecimal_inspect(VALUE self) return str; } +/* Returns self * 10**v without changing the precision. + * This method is currently for internal use. + * + * BigDecimal("0.123e10")._decimal_shift(20) #=> "0.123e30" + * BigDecimal("0.123e10")._decimal_shift(-20) #=> "0.123e-10" + */ +static VALUE +BigDecimal_decimal_shift(VALUE self, VALUE v) +{ + BDVALUE a, c; + ssize_t shift, exponentShift; + bool shiftDown; + size_t prec; + DECDIG ex, iex; + + a = GetBDValueMust(self); + shift = NUM2SSIZET(rb_to_int(v)); + + if (VpIsZero(a.real) || VpIsNaN(a.real) || VpIsInf(a.real) || shift == 0) return CheckGetValue(a); + + exponentShift = shift > 0 ? shift / BASE_FIG : (shift + 1) / BASE_FIG - 1; + shift -= exponentShift * BASE_FIG; + ex = 1; + for (int i = 0; i < shift; i++) ex *= 10; + shiftDown = a.real->frac[0] * (DECDIG_DBL)ex >= BASE; + iex = BASE / ex; + + prec = a.real->Prec + shiftDown; + c = NewZeroWrap(1, prec * BASE_FIG); + if (shift == 0) { + VpAsgn(c.real, a.real, 1); + } else if (shiftDown) { + DECDIG carry = 0; + exponentShift++; + for (size_t i = 0; i < a.real->Prec; i++) { + DECDIG v = a.real->frac[i]; + c.real->frac[i] = carry * ex + v / iex; + carry = v % iex; + } + c.real->frac[a.real->Prec] = carry * ex; + } else { + DECDIG carry = 0; + for (ssize_t i = a.real->Prec - 1; i >= 0; i--) { + DECDIG v = a.real->frac[i]; + c.real->frac[i] = v % iex * ex + carry; + carry = v / iex; + } + } + while (c.real->frac[prec - 1] == 0) prec--; + c.real->Prec = prec; + c.real->sign = a.real->sign; + c.real->exponent = a.real->exponent; + AddExponent(c.real, exponentShift); + RB_GC_GUARD(a.bigdecimal); + return CheckGetValue(c); +} + inline static int is_zero(VALUE x) { @@ -3564,6 +3622,7 @@ Init_bigdecimal(void) rb_define_method(rb_cBigDecimal, "infinite?", BigDecimal_IsInfinite, 0); rb_define_method(rb_cBigDecimal, "finite?", BigDecimal_IsFinite, 0); rb_define_method(rb_cBigDecimal, "truncate", BigDecimal_truncate, -1); + rb_define_method(rb_cBigDecimal, "_decimal_shift", BigDecimal_decimal_shift, 1); rb_define_method(rb_cBigDecimal, "_dump", BigDecimal_dump, -1); #ifdef BIGDECIMAL_USE_VP_TEST_METHODS @@ -3621,7 +3680,6 @@ enum op_sw { }; static int VpIsDefOP(Real *c, Real *a, Real *b, enum op_sw sw); -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); diff --git a/test/bigdecimal/test_bigdecimal.rb b/test/bigdecimal/test_bigdecimal.rb index 09f4e95f..77f5d278 100644 --- a/test/bigdecimal/test_bigdecimal.rb +++ b/test/bigdecimal/test_bigdecimal.rb @@ -1490,6 +1490,23 @@ def test_sqrt_minimum_precision assert_in_delta(BigDecimal('0.' + '0' * 25 + '3' * 100), x.sqrt(1), BigDecimal('1e-125')) end + def test_internal_use_decimal_shift + assert_positive_infinite_calculation { BigDecimal::INFINITY._decimal_shift(10) } + assert_negative_infinite_calculation { NEGATIVE_INFINITY._decimal_shift(10) } + assert_nan_calculation { BigDecimal::NAN._decimal_shift(10) } + assert_equal(BigDecimal(0), BigDecimal(0)._decimal_shift(10)) + [ + BigDecimal('-1234.56789'), + BigDecimal('123456789.012345678987654321'), + BigDecimal('123456789012345.678987654321') + ].each do |num| + (0..20).each do |shift| + assert_equal(num * 10**shift, num._decimal_shift(shift)) + assert_equal(num / 10**shift, num._decimal_shift(-shift)) + end + end + end + def test_fix x = BigDecimal("1.1") assert_equal(1, x.fix) From d6e2fffb48df7ec75264468377e53f27b8466549 Mon Sep 17 00:00:00 2001 From: tompng Date: Sat, 19 Jul 2025 18:10:56 +0900 Subject: [PATCH 2/2] Use _decimal_shift in sqrt, log and exp calculation --- lib/bigdecimal.rb | 55 ++++++++++++++++++++++------------------------- 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/lib/bigdecimal.rb b/lib/bigdecimal.rb index e1accd47..c8a00d77 100644 --- a/lib/bigdecimal.rb +++ b/lib/bigdecimal.rb @@ -197,7 +197,7 @@ def sqrt(prec) prec = [prec, n_digits].max ex = exponent / 2 - x = self.mult(BigDecimal("1e#{-ex * 2}"), n_significant_digits) + x = _decimal_shift(-2 * ex) y = BigDecimal(Math.sqrt(x.to_f)) precs = [prec + BigDecimal.double_fig] precs << 2 + precs.last / 2 while precs.last > BigDecimal.double_fig @@ -205,7 +205,7 @@ def sqrt(prec) y = y.add(x.div(y, p), p).div(2, p) end y = y.mult(1, limit) if limit - y.mult(BigDecimal("1e#{ex}"), precs.first) + y._decimal_shift(ex) end end @@ -240,7 +240,7 @@ def self.log(x, prec) if x > 10 || x < 0.1 log10 = log(BigDecimal(10), prec) exponent = x.exponent - x *= BigDecimal("1e#{-x.exponent}") + x = x._decimal_shift(-exponent) if x < 0.3 x *= 10 exponent -= 1 @@ -299,33 +299,30 @@ def self.exp(x, prec) return BigDecimal(1) if x.zero? return BigDecimal(1).div(exp(-x, prec), prec) if x < 0 - BigDecimal.save_limit do - BigDecimal.limit(0) - # exp(x * 10**cnt) = exp(x)**(10**cnt) - cnt = x > 1 ? x.exponent : 0 - prec2 = prec + BigDecimal.double_fig + cnt - x *= BigDecimal("1e-#{cnt}") - xn = BigDecimal(1) - y = BigDecimal(1) - - # Taylor series for exp(x) around 0 - 1.step do |i| - n = prec2 + xn.exponent - break if n <= 0 || xn.zero? - x = x.mult(1, n) - xn = xn.mult(x, n).div(i, n) - y = y.add(xn, prec2) - end - - # calculate exp(x * 10**cnt) from exp(x) - # exp(x * 10**k) = exp(x * 10**(k - 1)) ** 10 - cnt.times do - y2 = y.mult(y, prec2) - y5 = y2.mult(y2, prec2).mult(y, prec2) - y = y5.mult(y5, prec2) - end + # exp(x * 10**cnt) = exp(x)**(10**cnt) + cnt = x > 1 ? x.exponent : 0 + prec2 = prec + BigDecimal.double_fig + cnt + x = x._decimal_shift(-cnt) + xn = BigDecimal(1) + y = BigDecimal(1) + + # Taylor series for exp(x) around 0 + 1.step do |i| + n = prec2 + xn.exponent + break if n <= 0 || xn.zero? + x = x.mult(1, n) + xn = xn.mult(x, n).div(i, n) + y = y.add(xn, prec2) + end - y.mult(1, prec) + # calculate exp(x * 10**cnt) from exp(x) + # exp(x * 10**k) = exp(x * 10**(k - 1)) ** 10 + cnt.times do + y2 = y.mult(y, prec2) + y5 = y2.mult(y2, prec2).mult(y, prec2) + y = y5.mult(y5, prec2) end + + y.mult(1, prec) end end