Skip to content

Split internal extra calculation prec and BigDecimal.double_fig usage #486

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tompng merged 1 commit into from
Jan 30, 2026
+41 −39
Merged

Split internal extra calculation prec and BigDecimal.double_fig usage #486

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 7 additions & 4 deletions lib/bigdecimal.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,9 @@ def _decimal_shift(i) # :nodoc:

class BigDecimal
module Internal # :nodoc:
# Default extra precision for intermediate calculations
# This value is currently the same as BigDecimal.double_fig, but defined separately for future changes.
EXTRA_PREC = 16

# Coerce x to BigDecimal with the specified precision.
# TODO: some methods (example: BigMath.exp) require more precision than specified to coerce.
Expand Down Expand Up @@ -210,7 +213,7 @@ def power(y, prec = 0)
result_prec = prec.nonzero? || [x.n_significant_digits, y.n_significant_digits, BigDecimal.double_fig].max + BigDecimal.double_fig
result_prec = [result_prec, limit].min if prec.zero? && limit.nonzero?

prec2 = result_prec + BigDecimal.double_fig
prec2 = result_prec + BigDecimal::Internal::EXTRA_PREC

if y < 0
inv = x.power(-y, prec2)
Expand Down Expand Up @@ -266,7 +269,7 @@ def sqrt(prec)
ex = exponent / 2
x = _decimal_shift(-2 * ex)
y = BigDecimal(Math.sqrt(x.to_f), 0)
Internal.newton_loop(prec + BigDecimal.double_fig) do |p|
Internal.newton_loop(prec + BigDecimal::Internal::EXTRA_PREC) do |p|
y = y.add(x.div(y, p), p).div(2, p)
end
y._decimal_shift(ex).mult(1, prec)
Expand Down Expand Up @@ -301,7 +304,7 @@ def log(x, prec)
return BigDecimal::Internal.infinity_computation_result if x.infinite?
return BigDecimal(0) if x == 1

prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC

# Reduce x to near 1
if x > 1.01 || x < 0.99
Expand Down Expand Up @@ -352,7 +355,7 @@ def exp(x, prec)

# exp(x * 10**cnt) = exp(x)**(10**cnt)
cnt = x < -1 || x > 1 ? x.exponent : 0
prec2 = prec + BigDecimal.double_fig + cnt
prec2 = prec + BigDecimal::Internal::EXTRA_PREC + cnt
x = x._decimal_shift(-cnt)

# Decimal form of bit-burst algorithm
Expand Down
69 changes: 34 additions & 35 deletions lib/bigdecimal/math.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,19 +75,19 @@ def sqrt(x, prec)
private_class_method def _sin_periodic_reduction(x, prec, add_half_pi: false) # :nodoc:
return [1, x] if -Math::PI/2 <= x && x <= Math::PI/2 && !add_half_pi

mod_prec = prec + BigDecimal.double_fig
pi_extra_prec = [x.exponent, 0].max + BigDecimal.double_fig
mod_prec = prec + BigDecimal::Internal::EXTRA_PREC
pi_extra_prec = [x.exponent, 0].max + BigDecimal::Internal::EXTRA_PREC
while true
pi = PI(mod_prec + pi_extra_prec)
half_pi = pi / 2
div, mod = (add_half_pi ? x + pi : x + half_pi).divmod(pi)
mod -= half_pi
if mod.zero? || mod_prec + mod.exponent <= 0
# mod is too small to estimate required pi precision
mod_prec = mod_prec * 3 / 2 + BigDecimal.double_fig
mod_prec = mod_prec * 3 / 2 + BigDecimal::Internal::EXTRA_PREC
elsif mod_prec + mod.exponent < prec
# Estimate required precision of pi
mod_prec = prec - mod.exponent + BigDecimal.double_fig
mod_prec = prec - mod.exponent + BigDecimal::Internal::EXTRA_PREC
else
return [div % 2 == 0 ? 1 : -1, mod.mult(1, prec)]
end
Expand Down Expand Up @@ -145,9 +145,7 @@ def cbrt(x, prec)
ex = x.exponent / 3
x = x._decimal_shift(-3 * ex)
y = BigDecimal(Math.cbrt(x.to_f), 0)
precs = [prec + BigDecimal.double_fig]
precs << 2 + precs.last / 2 while precs.last > BigDecimal.double_fig
precs.reverse_each do |p|
BigDecimal::Internal.newton_loop(prec + BigDecimal::Internal::EXTRA_PREC) do |p|
y = (2 * y + x.div(y, p).div(y, p)).div(3, p)
end
y._decimal_shift(ex).mult(neg ? -1 : 1, prec)
Expand All @@ -168,7 +166,7 @@ def hypot(x, y, prec)
y = BigDecimal::Internal.coerce_to_bigdecimal(y, prec, :hypot)
return BigDecimal::Internal.nan_computation_result if x.nan? || y.nan?
return BigDecimal::Internal.infinity_computation_result if x.infinite? || y.infinite?
prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
sqrt(x.mult(x, prec2) + y.mult(y, prec2), prec)
end

Expand All @@ -187,7 +185,7 @@ def sin(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :sin)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sin)
return BigDecimal::Internal.nan_computation_result if x.infinite? || x.nan?
n = prec + BigDecimal.double_fig
n = prec + BigDecimal::Internal::EXTRA_PREC
sign, x = _sin_periodic_reduction(x, n)
_sin_around_zero(x, n).mult(sign, prec)
end
Expand All @@ -207,7 +205,7 @@ def cos(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :cos)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cos)
return BigDecimal::Internal.nan_computation_result if x.infinite? || x.nan?
n = prec + BigDecimal.double_fig
n = prec + BigDecimal::Internal::EXTRA_PREC
sign, x = _sin_periodic_reduction(x, n, add_half_pi: true)
_sin_around_zero(x, n).mult(sign, prec)
end
Expand All @@ -228,7 +226,8 @@ def cos(x, prec)
#
def tan(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :tan)
sin(x, prec + BigDecimal.double_fig).div(cos(x, prec + BigDecimal.double_fig), prec)
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
sin(x, prec2).div(cos(x, prec2), prec)
end

# call-seq:
Expand All @@ -248,7 +247,7 @@ def asin(x, prec)
raise Math::DomainError, "Out of domain argument for asin" if x < -1 || x > 1
return BigDecimal::Internal.nan_computation_result if x.nan?

prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
cos = (1 - x**2).sqrt(prec2)
if cos.zero?
PI(prec2).div(x > 0 ? 2 : -2, prec)
Expand All @@ -274,7 +273,7 @@ def acos(x, prec)
raise Math::DomainError, "Out of domain argument for acos" if x < -1 || x > 1
return BigDecimal::Internal.nan_computation_result if x.nan?

prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
return (PI(prec2) / 2).sub(asin(x, prec2), prec) if x < 0
return PI(prec2).div(2, prec) if x.zero?

Expand All @@ -297,7 +296,7 @@ def atan(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :atan)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atan)
return BigDecimal::Internal.nan_computation_result if x.nan?
n = prec + BigDecimal.double_fig
n = prec + BigDecimal::Internal::EXTRA_PREC
return PI(n).div(2 * x.infinite?, prec) if x.infinite?

x = -x if neg = x < 0
Expand Down Expand Up @@ -341,7 +340,7 @@ def atan2(y, x, prec)

y = -y if neg = y < 0
xlarge = y.abs < x.abs
prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
if x > 0
v = xlarge ? atan(y.div(x, prec2), prec) : PI(prec2) / 2 - atan(x.div(y, prec2), prec2)
else
Expand All @@ -367,7 +366,7 @@ def sinh(x, prec)
return BigDecimal::Internal.nan_computation_result if x.nan?
return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite?

prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
prec2 -= x.exponent if x.exponent < 0
e = exp(x, prec2)
(e - BigDecimal(1).div(e, prec2)).div(2, prec)
Expand All @@ -390,7 +389,7 @@ def cosh(x, prec)
return BigDecimal::Internal.nan_computation_result if x.nan?
return BigDecimal::Internal.infinity_computation_result if x.infinite?

prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
e = exp(x, prec2)
(e + BigDecimal(1).div(e, prec2)).div(2, prec)
end
Expand All @@ -412,7 +411,7 @@ def tanh(x, prec)
return BigDecimal::Internal.nan_computation_result if x.nan?
return BigDecimal(x.infinite?) if x.infinite?

prec2 = prec + BigDecimal.double_fig + [-x.exponent, 0].max
prec2 = prec + BigDecimal::Internal::EXTRA_PREC + [-x.exponent, 0].max
e = exp(x, prec2)
einv = BigDecimal(1).div(e, prec2)
(e - einv).div(e + einv, prec)
Expand All @@ -436,7 +435,7 @@ def asinh(x, prec)
return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite?
return -asinh(-x, prec) if x < 0

sqrt_prec = prec + [-x.exponent, 0].max + BigDecimal.double_fig
sqrt_prec = prec + [-x.exponent, 0].max + BigDecimal::Internal::EXTRA_PREC
log(x + sqrt(x**2 + 1, sqrt_prec), prec)
end

Expand All @@ -458,7 +457,7 @@ def acosh(x, prec)
return BigDecimal::Internal.infinity_computation_result if x.infinite?
return BigDecimal::Internal.nan_computation_result if x.nan?

log(x + sqrt(x**2 - 1, prec + BigDecimal.double_fig), prec)
log(x + sqrt(x**2 - 1, prec + BigDecimal::Internal::EXTRA_PREC), prec)
end

# call-seq:
Expand All @@ -480,7 +479,7 @@ def atanh(x, prec)
return BigDecimal::Internal.infinity_computation_result if x == 1
return -BigDecimal::Internal.infinity_computation_result if x == -1

prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
(log(x + 1, prec2) - log(1 - x, prec2)).div(2, prec)
end

Expand All @@ -505,10 +504,10 @@ def log2(x, prec)
return BigDecimal::Internal.nan_computation_result if x.nan?
return BigDecimal::Internal.infinity_computation_result if x.infinite? == 1

prec2 = prec + BigDecimal.double_fig * 3 / 2
prec2 = prec + BigDecimal::Internal::EXTRA_PREC * 3 / 2
v = log(x, prec2).div(log(BigDecimal(2), prec2), prec2)
# Perform half-up rounding to calculate log2(2**n)==n correctly in every rounding mode
v = v.round(prec + BigDecimal.double_fig - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP)
v = v.round(prec + BigDecimal::Internal::EXTRA_PREC - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP)
v.mult(1, prec)
end

Expand All @@ -533,10 +532,10 @@ def log10(x, prec)
return BigDecimal::Internal.nan_computation_result if x.nan?
return BigDecimal::Internal.infinity_computation_result if x.infinite? == 1

prec2 = prec + BigDecimal.double_fig * 3 / 2
prec2 = prec + BigDecimal::Internal::EXTRA_PREC * 3 / 2
v = log(x, prec2).div(log(BigDecimal(10), prec2), prec2)
# Perform half-up rounding to calculate log10(10**n)==n correctly in every rounding mode
v = v.round(prec + BigDecimal.double_fig - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP)
v = v.round(prec + BigDecimal::Internal::EXTRA_PREC - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP)
v.mult(1, prec)
end

Expand Down Expand Up @@ -573,9 +572,9 @@ def expm1(x, prec)
if x < -1
# log10(exp(x)) = x * log10(e)
lg_e = 0.4342944819032518
exp_prec = prec + (lg_e * x).ceil + BigDecimal.double_fig
exp_prec = prec + (lg_e * x).ceil + BigDecimal::Internal::EXTRA_PREC
elsif x < 1
exp_prec = prec - x.exponent + BigDecimal.double_fig
exp_prec = prec - x.exponent + BigDecimal::Internal::EXTRA_PREC
else
exp_prec = prec
end
Expand Down Expand Up @@ -613,7 +612,7 @@ def erf(x, prec)
return BigDecimal(1).sub(erfc, prec) if erfc
end

prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
x_smallprec = x.mult(1, Integer.sqrt(prec2) / 2)
# Taylor series of x with small precision is fast
erf1 = _erf_taylor(x_smallprec, BigDecimal(0), BigDecimal(0), prec2)
Expand All @@ -637,7 +636,7 @@ def erfc(x, prec)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erfc)
return BigDecimal::Internal.nan_computation_result if x.nan?
return BigDecimal(1 - x.infinite?) if x.infinite?
return BigDecimal(1).sub(erf(x, prec + BigDecimal.double_fig), prec) if x < 0
return BigDecimal(1).sub(erf(x, prec + BigDecimal::Internal::EXTRA_PREC), prec) if x < 0
return BigDecimal(0) if x > 5000000000 # erfc(5000000000) < 1e-10000000000000000000 (underflow)

if x >= 8
Expand All @@ -648,7 +647,7 @@ def erfc(x, prec)
# erfc(x) = 1 - erf(x) < exp(-x**2)/x/sqrt(pi)
# Precision of erf(x) needs about log10(exp(-x**2)) extra digits
log10 = 2.302585092994046
high_prec = prec + BigDecimal.double_fig + (x.ceil**2 / log10).ceil
high_prec = prec + BigDecimal::Internal::EXTRA_PREC + (x.ceil**2 / log10).ceil
BigDecimal(1).sub(erf(x, high_prec), prec)
end

Expand Down Expand Up @@ -697,7 +696,7 @@ def erfc(x, prec)
# Using Stirling's approximation, we can simplify this condition to:
# sqrt(2)/2 + k*log(k) - k - 2*k*log(x) < -prec*log(10)
# and the left side is minimized when k = x**2.
prec += BigDecimal.double_fig
prec += BigDecimal::Internal::EXTRA_PREC
xf = x.to_f
kmax = (1..(xf ** 2).floor).bsearch do |k|
Math.log(2) / 2 + k * Math.log(k) - k - 2 * k * Math.log(xf) < -prec * Math.log(10)
Expand Down Expand Up @@ -726,7 +725,7 @@ def erfc(x, prec)
def gamma(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma)
prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
if x < 0.5
raise Math::DomainError, 'Numerical argument is out of domain - gamma' if x.frac.zero?

Expand Down Expand Up @@ -754,7 +753,7 @@ def gamma(x, prec)
def lgamma(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma)
prec2 = prec + BigDecimal.double_fig
prec2 = prec + BigDecimal::Internal::EXTRA_PREC
if x < 0.5
return [BigDecimal::INFINITY, 1] if x.frac.zero?

Expand All @@ -763,7 +762,7 @@ def lgamma(x, prec)
pi = PI(prec2)
sin = _sinpix(x, pi, prec2)
log_gamma = BigMath.log(pi, prec2).sub(lgamma(1 - x, prec2).first + BigMath.log(sin.abs, prec2), prec)
return [log_gamma, sin > 0 ? 1 : -1] if prec2 + log_gamma.exponent > prec + BigDecimal.double_fig
return [log_gamma, sin > 0 ? 1 : -1] if prec2 + log_gamma.exponent > prec + BigDecimal::Internal::EXTRA_PREC

# Retry with higher precision if loss of significance is too large
prec2 = prec2 * 3 / 2
Expand Down Expand Up @@ -896,7 +895,7 @@ def ldexp(x, exponent)
def PI(prec)
# Gauss–Legendre algorithm
prec = BigDecimal::Internal.coerce_validate_prec(prec, :PI)
n = prec + BigDecimal.double_fig
n = prec + BigDecimal::Internal::EXTRA_PREC
a = BigDecimal(1)
b = BigDecimal(0.5, 0).sqrt(n)
s = BigDecimal(0.25, 0)
Expand Down
Loading