Skip to content

Implement BigMath.erf(x, prec) and BigMath.erfc(x, prec) #357

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
mrkn merged 1 commit into from
Dec 10, 2025
Merged

Implement BigMath.erf(x, prec) and BigMath.erfc(x, prec) #357

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
130 changes: 130 additions & 0 deletions lib/bigdecimal/math.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@
# log10(x, prec)
# log1p(x, prec)
# expm1(x, prec)
# erf (x, prec)
# erfc(x, prec)
# PI (prec)
# E (prec) == exp(1.0,prec)
#
Expand Down Expand Up @@ -568,6 +570,134 @@ def expm1(x, prec)
exp_prec > 0 ? exp(x, exp_prec).sub(1, prec) : BigDecimal(-1)
end

# call-seq:
# erf(decimal, numeric) -> BigDecimal
#
# Computes the error function of +decimal+ to the specified number of digits of
# precision, +numeric+.
#
# If +decimal+ is NaN, returns NaN.
#
# BigMath.erf(BigDecimal('1'), 32).to_s
# #=> "0.84270079294971486934122063508261e0"
#
def erf(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :erf)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erf)
return BigDecimal::Internal.nan_computation_result if x.nan?
return BigDecimal(x.infinite?) if x.infinite?
return BigDecimal(0) if x == 0
return -erf(-x, prec) if x < 0
return BigDecimal(1) if x > 5000000000 # erf(5000000000) > 1 - 1e-10000000000000000000

if x > 8
xf = x.to_f
log10_erfc = -xf ** 2 / Math.log(10) - Math.log10(xf * Math::PI ** 0.5)
erfc_prec = [prec + log10_erfc.ceil, 1].max
erfc = _erfc_asymptotic(x, erfc_prec)
return BigDecimal(1).sub(erfc, prec) if erfc
end

prec2 = prec + BigDecimal.double_fig
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)
# Taylor series converges quickly for small x
_erf_taylor(x - x_smallprec, x_smallprec, erf1, prec2).mult(1, prec)
end

# call-seq:
# erfc(decimal, numeric) -> BigDecimal
#
# Computes the complementary error function of +decimal+ to the specified number of digits of
# precision, +numeric+.
#
# If +decimal+ is NaN, returns NaN.
#
# BigMath.erfc(BigDecimal('10'), 32).to_s
# #=> "0.20884875837625447570007862949578e-44"
#
def erfc(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :erfc)
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), prec) if x < 0
return BigDecimal(0) if x > 5000000000 # erfc(5000000000) < 1e-10000000000000000000 (underflow)

if x >= 8
y = _erfc_asymptotic(x, prec)
return y.mult(1, prec) if y
end

# 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
BigDecimal(1).sub(erf(x, high_prec), prec)
end

# Calculates erf(x + a)
private_class_method def _erf_taylor(x, a, erf_a, prec) # :nodoc:
return erf_a if x.zero?
# Let f(x+a) = erf(x+a)*exp((x+a)**2)*sqrt(pi)/2
# = c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ...
# f'(x+a) = 1+2*(x+a)*f(x+a)
# f'(x+a) = c1 + 2*c2*x + 3*c3*x**2 + 4*c4*x**3 + 5*c5*x**4 + ...
# = 1+2*(x+a)*(c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ...)
# therefore,
# c0 = f(a)
# c1 = 2 * a * c0 + 1
# c2 = (2 * c0 + 2 * a * c1) / 2
# c3 = (2 * c1 + 2 * a * c2) / 3
# c4 = (2 * c2 + 2 * a * c3) / 4
#
# All coefficients are positive when a >= 0

scale = BigDecimal(2).div(sqrt(PI(prec), prec), prec)
c_prev = erf_a.div(scale.mult(exp(-a*a, prec), prec), prec)
c_next = (2 * a * c_prev).add(1, prec).mult(x, prec)
sum = c_prev.add(c_next, prec)

2.step do |k|
c = (c_prev.mult(x, prec) + a * c_next).mult(2, prec).mult(x, prec).div(k, prec)
sum = sum.add(c, prec)
c_prev, c_next = c_next, c
break if [c_prev, c_next].all? { |c| c.zero? || (c.exponent < sum.exponent - prec) }
end
value = sum.mult(scale.mult(exp(-(x + a).mult(x + a, prec), prec), prec), prec)
value > 1 ? BigDecimal(1) : value
end

private_class_method def _erfc_asymptotic(x, prec) # :nodoc:
# Let f(x) = erfc(x)*sqrt(pi)*exp(x**2)/2
# f(x) satisfies the following differential equation:
# 2*x*f(x) = f'(x) + 1
# From the above equation, we can derive the following asymptotic expansion:
# f(x) = (0..kmax).sum { (-1)**k * (2*k)! / 4**k / k! / x**(2*k)) } / x

# This asymptotic expansion does not converge.
# But if there is a k that satisfies (2*k)! / 4**k / k! / x**(2*k) < 10**(-prec),
# It is enough to calculate erfc within the given precision.
# 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
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)
end
return unless kmax

sum = BigDecimal(1)
x2 = x.mult(x, prec)
d = BigDecimal(1)
(1..kmax).each do |k|
d = d.div(x2, prec).mult(1 - 2 * k, prec).div(2, prec)
sum = sum.add(d, prec)
end
sum.div(exp(x2, prec).mult(PI(prec).sqrt(prec), prec), prec).div(x, prec)
end

# call-seq:
# PI(numeric) -> BigDecimal
Expand Down
58 changes: 58 additions & 0 deletions test/bigdecimal/test_bigmath.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,8 @@ def test_consistent_precision_acceptance
assert_consistent_precision_acceptance {|prec| BigMath.log10(x, prec) }
assert_consistent_precision_acceptance {|prec| BigMath.log1p(x, prec) }
assert_consistent_precision_acceptance {|prec| BigMath.expm1(x, prec) }
assert_consistent_precision_acceptance {|prec| BigMath.erf(x, prec) }
assert_consistent_precision_acceptance {|prec| BigMath.erfc(x, prec) }

assert_consistent_precision_acceptance {|prec| BigMath.E(prec) }
assert_consistent_precision_acceptance {|prec| BigMath.PI(prec) }
Expand Down Expand Up @@ -112,6 +114,8 @@ def test_coerce_argument
assert_equal(log10(bd, N), log10(f, N))
assert_equal(log1p(bd, N), log1p(f, N))
assert_equal(expm1(bd, N), expm1(f, N))
assert_equal(erf(bd, N), erf(f, N))
assert_equal(erfc(bd, N), erfc(f, N))
end

def test_sqrt
Expand Down Expand Up @@ -469,4 +473,58 @@ def test_expm1
assert_in_exact_precision(exp(BigDecimal("1.23e-10"), 120) - 1, expm1(BigDecimal("1.23e-10"), 100), 100)
assert_in_exact_precision(exp(123, 120) - 1, expm1(BigDecimal("123"), 100), 100)
end

def test_erf
[-0.5, 0.1, 0.3, 2.1, 3.3].each do |x|
assert_in_epsilon(Math.erf(x), BigMath.erf(BigDecimal(x.to_s), N))
end
assert_equal(1, BigMath.erf(PINF, N))
assert_equal(-1, BigMath.erf(MINF, N))
assert_equal(1, BigMath.erf(BigDecimal(1000), 100))
assert_equal(-1, BigMath.erf(BigDecimal(-1000), 100))
assert_not_equal(1, BigMath.erf(BigDecimal(10), 45))
assert_not_equal(1, BigMath.erf(BigDecimal(15), 100))
assert_equal(1, BigMath.erf(BigDecimal('1e400'), 10))
assert_equal(-1, BigMath.erf(BigDecimal('-1e400'), 10))
assert_equal(BigMath.erf(BigDecimal('1e-300'), N) * BigDecimal('1e-100'), BigMath.erf(BigDecimal('1e-400'), N))
assert_equal(
BigDecimal("0.9953222650189527341620692563672529286108917970400600767383523262004372807199951773676290080196806805"),
BigMath.erf(BigDecimal("2"), 100)
)
assert_converge_in_precision {|n| BigMath.erf(BigDecimal("1e-30"), n) }
assert_converge_in_precision {|n| BigMath.erf(BigDecimal("0.3"), n) }
assert_converge_in_precision {|n| BigMath.erf(SQRT2, n) }
end

def test_erfc
[-0.5, 0.1, 0.3, 2.1, 3.3].each do |x|
assert_in_epsilon(Math.erfc(x), BigMath.erfc(BigDecimal(x.to_s), N))
end
assert_equal(0, BigMath.erfc(PINF, N))
assert_equal(2, BigMath.erfc(MINF, N))
assert_equal(0, BigMath.erfc(BigDecimal('1e400'), 10))
assert_equal(2, BigMath.erfc(BigDecimal('-1e400'), 10))
assert_equal(1, BigMath.erfc(BigDecimal('1e-400'), N))

# erfc with taylor series
assert_equal(
BigDecimal("2.088487583762544757000786294957788611560818119321163727012213713938174695833440290610766384285723554e-45"),
BigMath.erfc(BigDecimal("10"), 100)
)
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(0.3), n) }
assert_converge_in_precision {|n| BigMath.erfc(SQRT2, n) }
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(8), n) }
# erfc with asymptotic expansion
assert_equal(
BigDecimal("1.896961059966276509268278259713415434936907563929186183462834752900411805205111886605256690776760041e-697"),
BigMath.erfc(BigDecimal("40"), 100)
)
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(30), n) }
assert_converge_in_precision {|n| BigMath.erfc(30 * SQRT2, n) }
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(50), n) }
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(60000), n) }
# Near crossover point between taylor series and asymptotic expansion around prec=150
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(19.5), n) }
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(20.5), n) }
end
end
2 changes: 2 additions & 0 deletions test/bigdecimal/test_jruby.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,8 @@ def test_bigmath
assert_in_delta(Math.log10(3), BigMath.log10(BigDecimal(3), N))
assert_in_delta(Math.log1p(0.1), BigMath.log1p(BigDecimal('0.1'), N)) if defined? Math.log1p
assert_in_delta(Math.expm1(0.1), BigMath.expm1(BigDecimal('0.1'), N)) if defined? Math.expm1
assert_in_delta(Math.erf(1), BigMath.erf(BigDecimal(1), N))
assert_in_delta(Math.erfc(10), BigMath.erfc(BigDecimal(10), N))
assert_in_delta(Math::PI, BigMath.PI(N))
assert_in_delta(Math::E, BigMath.E(N))
end
Expand Down