Implement BigMath.cbrt

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,34 @@ 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
+
+    ex = x.exponent / 3
+    x *= BigDecimal("1e#{-ex * 3}")
+    y = BigDecimal(Math.cbrt(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 = (2 * y + x.div(y, p).div(y, p)).div(3, p)
+    end
+    y * BigDecimal("1e#{ex}")
+  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"))