Implement BigMath.cbrt

Implement BigMath version of Math.cbrt(x).
Added internal-use method BigMath.int_cbrt(n): cube root version of Integer.sqrt(n)

Author: tompng

lib/bigdecimal/math.rb

@@ -5,6 +5,7 @@
 #--
 # Contents:
 #   sqrt(x, prec)
+#   cbrt(x, prec)
 #   hypot(x, y, prec)
 #   sin (x, prec)
 #   cos (x, prec)
@@ -46,6 +47,61 @@ def sqrt(x, prec)
     x.sqrt(prec)
   end
 
+  # call-seq:
+  #   cbrt(decimal, numeric) -> BigDecimal
+  #
+  # Computes the cube root of +decimal+ to the specified number of digits of
+  # precision, +numeric+.
+  #
+  #   BigMath.cbrt(BigDecimal('2'), 16).to_s
+  #   #=> "0.125992104989487316476721060727822e1"
+  #
+  def cbrt(x, prec)
+    raise ArgumentError, "Zero or negative precision for cbrt" if prec <= 0
+    return x if x.zero? || x.infinite? || x.nan?
+    return -cbrt(-x, prec) if x < 0
+
+    n_digits = x.n_significant_digits
+    prec = [prec, n_digits].max
+
+    if n_digits < prec / 2
+      # Fast path for cbrt(8e150) => 2e50
+      ex = (n_digits - x.exponent + 2) / 3
+      n = (x * BigDecimal("1e#{3 * ex}")).to_i
+      cbrt = _int_cbrt(n)
+      return BigDecimal(cbrt) * BigDecimal("1e#{-ex}") if cbrt**3 == n
+    end
+
+    ex = prec + BigDecimal.double_fig - x.exponent / 3
+    cbrt = _int_cbrt(x * BigDecimal("1e#{3 * ex}"))
+    BigDecimal(cbrt) * BigDecimal("1e#{-ex}")
+  end
+
+  # Private method used internally by `cbrt`.
+  # Cube root version of `Intger.sqrt(n)`
+  # Returns the largest integer whose cube is less than or equal to n if n is positive.
+  private_class_method def _int_cbrt(n)
+    n = n.to_i
+    return -_int_cbrt(-n) if n < 0
+
+    if n <= 0xffffffff
+      v = Math.cbrt(n).floor
+    else
+      shift = (n.bit_length - 1) / 6
+      v = _int_cbrt(n >> (3 * shift)) << shift
+    end
+
+    v = (2 * v + n / v / v) / 3
+    v2 = v * v
+    v3 = v2 * v
+    while v3 > n
+      v3 -= 3 * v2 - 3 * v + 1
+      v2 -= 2 * v - 1
+      v -= 1
+    end
+    v
+  end
+
   # call-seq:
   #   hypot(x, y, numeric) -> BigDecimal
   #

test/bigdecimal/test_bigmath.rb

@@ -34,6 +34,26 @@ def test_sqrt
     assert_relative_precision {|n| sqrt(BigDecimal("2e50"), n) }
   end
 
+  def test_cbrt
+    assert_equal(1234, cbrt(BigDecimal(1234**3), N))
+    assert_equal(-12345, cbrt(BigDecimal(-12345**3), N))
+    assert_equal(123 ** 456, cbrt(BigDecimal(123 ** 456) ** 3, N))
+    assert_equal(0, cbrt(BigDecimal("0"), N))
+    assert_equal(0, cbrt(BigDecimal("-0"), N))
+    assert_equal(PINF, cbrt(PINF, N))
+    assert_equal(MINF, cbrt(MINF, N))
+
+    assert_in_delta(SQRT2, cbrt(SQRT2 ** 3, 100), BigDecimal("1e-100"))
+    assert_in_delta(SQRT3, cbrt(SQRT3 ** 3, 100), BigDecimal("1e-100"))
+    assert_equal(BigDecimal("3e50"), cbrt(BigDecimal("27e150"), N))
+    assert_equal(BigDecimal("-4e50"), cbrt(BigDecimal("-64e150"), N))
+    assert_in_epsilon(Math.cbrt(28e150), cbrt(BigDecimal("28e150"), N))
+    assert_in_epsilon(Math.cbrt(27e151), cbrt(BigDecimal("27e151"), N))
+    assert_relative_precision {|n| cbrt(BigDecimal("2"), n) }
+    assert_relative_precision {|n| cbrt(BigDecimal("2e-50"), n) }
+    assert_relative_precision {|n| cbrt(BigDecimal("2e50"), n) }
+  end
+
   def test_hypot
     assert_in_delta(SQRT2, hypot(BigDecimal("1"), BigDecimal("1"), 100), BigDecimal("1e-100"))
     assert_in_delta(SQRT5, hypot(SQRT2, SQRT3, 100), BigDecimal("1e-100"))