Implement BigDecimal#_decimal_shift for internal use (#324)

tompng · web-flow · commit b946e79a3d45 · 2025-09-02T23:57:28.000+09:00
* Implement BigDecimal#_decimal_shift for internal use

* Use _decimal_shift in sqrt, log and exp calculation
diff --git 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/lib/bigdecimal.rb b/lib/bigdecimal.rb
@@ -197,15 +197,15 @@ 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
     precs.reverse_each do |p|
       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
diff --git 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)
