Skip to content

Commit 09dce1c

Browse files
committed
Improve performance of x**y when y is a huge value
When y.exponent is several thousand or more, x**y was slow because repeated squaring algorithm requires several thousands of multiplications. Skip repeated squaring algorithm in such case. Needed to calaculate (1+1/n).power(n, prec)
1 parent cb2458b commit 09dce1c

File tree

2 files changed

+28
-11
lines changed

2 files changed

+28
-11
lines changed

lib/bigdecimal.rb

Lines changed: 23 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -165,20 +165,32 @@ def power(y, prec = nil)
165165
return BigDecimal(1).div(inv, prec)
166166
end
167167

168-
int_part = y.fix.to_i
168+
if y.exponent > Math.log(prec) * 5 + 20
169+
# x**int_part calculation is slow when y.exponent is large, so skip it
170+
int_part = 0
171+
rest_part = y
172+
else
173+
int_part = y.fix.to_i
174+
rest_part = frac_part
175+
end
176+
169177
prec2 = prec + BigDecimal.double_fig
170-
pow_prec = prec2 + (int_part > 0 ? y.exponent : 0)
171178
ans = BigDecimal(1)
172-
n = 1
173-
xn = x
174-
while true
175-
ans = ans.mult(xn, pow_prec) if int_part.allbits?(n)
176-
n <<= 1
177-
break if n > int_part
178-
xn = xn.mult(xn, pow_prec)
179+
180+
if int_part > 0
181+
pow_prec = prec2 + y.exponent
182+
n = 1
183+
xn = x
184+
while true
185+
ans = ans.mult(xn, pow_prec) if int_part.allbits?(n)
186+
n <<= 1
187+
break if n > int_part
188+
xn = xn.mult(xn, pow_prec)
189+
end
179190
end
180-
unless frac_part.zero?
181-
ans = ans.mult(BigMath.exp(BigMath.log(x, prec2).mult(frac_part, prec2), prec2), prec2)
191+
192+
unless rest_part.zero?
193+
ans = ans.mult(BigMath.exp(BigMath.log(x, prec2).mult(rest_part, prec2), prec2), prec2)
182194
end
183195
ans.mult(1, prec)
184196
end

test/bigdecimal/test_bigdecimal.rb

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1991,6 +1991,11 @@ def test_power_with_rational
19911991
assert_in_epsilon(z2, x2 ** y, 1e-99)
19921992
end
19931993

1994+
def test_power_with_huge_value
1995+
n = BigDecimal('7e+10000')
1996+
assert_equal(BigMath.exp(1, 100), (1 + BigDecimal(1).div(n, 120)).power(n, 100))
1997+
end
1998+
19941999
def test_power_precision
19952000
x = BigDecimal("1.41421356237309504880168872420969807856967187537695")
19962001
y = BigDecimal("3.14159265358979323846264338327950288419716939937511")

0 commit comments

Comments
 (0)