Fix precision of x.power(y, prec) when the result is nearly infinity (#439)

Author: tompng

lib/bigdecimal.rb

@@ -182,6 +182,13 @@ def power(y, prec = nil)
       end
       ans.mult(1, prec)
     else
+      if x > 1
+        # To calculate exp(z, prec), z needs prec+max(z.exponent, 0) precision if z > 0.
+        # Estimate (y*log(x)).exponent
+        logx_exponent = x < 2 ? (x - 1).exponent : Math.log10(x.exponent).round
+        ylogx_exponent = y.exponent + logx_exponent
+        prec2 += [ylogx_exponent, 0].max
+      end
       BigMath.exp(BigMath.log(x, prec2).mult(y, prec2), prec)
     end
   end

test/bigdecimal/test_bigdecimal.rb

@@ -1996,6 +1996,22 @@ def test_power_with_huge_value
     assert_equal(BigMath.exp(1, 100), (1 + BigDecimal(1).div(n, 120)).power(n, 100))
   end
 
+  def test_power_almost_infinite
+    omit if EXPONENT_MAX < 9e+18
+
+    x = BigDecimal(1.9)
+    y = BigDecimal('1e+18') + 0.1
+    assert_equal(x.power(y, 120).mult(1, 100), x.power(y, 100))
+
+    x = BigDecimal("9e+#{10**17}").div(7, 100)
+    y = BigDecimal(1.9)
+    assert_equal(x.power(y, 120).mult(1, 100), x.power(y, 100))
+
+    x = BigDecimal("9e+#{10**9}").div(7, 100)
+    y = BigDecimal('1e+8') + 0.1
+    assert_equal(x.power(y, 120).mult(1, 100), x.power(y, 100))
+  end
+
   def test_power_precision
     x = BigDecimal("1.41421356237309504880168872420969807856967187537695")
     y = BigDecimal("3.14159265358979323846264338327950288419716939937511")