Skip to content

Commit bb3cd2f

Browse files
tompngtompng
committed
Implement BigMath.erf(x, prec) and BigMath.erfc(x, prec)
Uses asymptotic expansion of erfc if possible and fallback to taylor series of erf
1 parent 150d01d commit bb3cd2f

File tree

3 files changed

+190
-0
lines changed

3 files changed

+190
-0
lines changed

‎lib/bigdecimal/math.rb‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@
2424
# log10(x, prec)
2525
# log1p(x, prec)
2626
# expm1(x, prec)
27+
# erf (x, prec)
28+
# erfc(x, prec)
2729
# PI (prec)
2830
# E (prec) == exp(1.0,prec)
2931
#
@@ -568,6 +570,134 @@ def expm1(x, prec)
568570
exp_prec > 0 ? exp(x, exp_prec).sub(1, prec) : BigDecimal(-1)
569571
end
570572

573+
# call-seq:
574+
# erf(decimal, numeric) -> BigDecimal
575+
#
576+
# Computes the error function of +decimal+ to the specified number of digits of
577+
# precision, +numeric+.
578+
#
579+
# If +decimal+ is NaN, returns NaN.
580+
#
581+
# BigMath.erf(BigDecimal('1'), 32).to_s
582+
# #=> "0.84270079294971486934122063508261e0"
583+
#
584+
def erf(x, prec)
585+
prec = BigDecimal::Internal.coerce_validate_prec(prec, :erf)
586+
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erf)
587+
return BigDecimal::Internal.nan_computation_result if x.nan?
588+
return BigDecimal(x.infinite?) if x.infinite?
589+
return BigDecimal(0) if x == 0
590+
return -erf(-x, prec) if x < 0
591+
return BigDecimal(1) if x > 5000000000 # erf(5000000000) > 1 - 1e-10000000000000000000
592+
593+
if x > 8
594+
xf = x.to_f
595+
log10_erfc = -xf ** 2 / Math.log(10) - Math.log10(xf * Math::PI ** 0.5)
596+
erfc_prec = [prec + log10_erfc.ceil, 1].max
597+
erfc = _erfc_asymptotic(x, erfc_prec)
598+
return BigDecimal(1).sub(erfc, prec) if erfc
599+
end
600+
601+
prec2 = prec + BigDecimal.double_fig
602+
x_smallprec = x.mult(1, Integer.sqrt(prec2) / 2)
603+
# Taylor series of x with small precision is fast
604+
erf1 = _erf_taylor(x_smallprec, BigDecimal(0), BigDecimal(0), prec2)
605+
# Taylor series converges quickly for small x
606+
_erf_taylor(x - x_smallprec, x_smallprec, erf1, prec2).mult(1, prec)
607+
end
608+
609+
# call-seq:
610+
# erfc(decimal, numeric) -> BigDecimal
611+
#
612+
# Computes the complementary error function of +decimal+ to the specified number of digits of
613+
# precision, +numeric+.
614+
#
615+
# If +decimal+ is NaN, returns NaN.
616+
#
617+
# BigMath.erfc(BigDecimal('10'), 32).to_s
618+
# #=> "0.20884875837625447570007862949578e-44"
619+
#
620+
def erfc(x, prec)
621+
prec = BigDecimal::Internal.coerce_validate_prec(prec, :erfc)
622+
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erfc)
623+
return BigDecimal::Internal.nan_computation_result if x.nan?
624+
return BigDecimal(1 - x.infinite?) if x.infinite?
625+
return BigDecimal(1).sub(erf(x, prec), prec) if x < 0
626+
return BigDecimal(0) if x > 5000000000 # erfc(5000000000) < 1e-10000000000000000000 (underflow)
627+
628+
if x >= 8
629+
y = _erfc_asymptotic(x, prec)
630+
return y.mult(1, prec) if y
631+
end
632+
633+
# erfc(x) = 1 - erf(x) < exp(-x**2)/x/sqrt(pi)
634+
# Precision of erf(x) needs about log10(exp(-x**2)) extra digits
635+
log10 = 2.302585092994046
636+
high_prec = prec + BigDecimal.double_fig + (x.ceil**2 / log10).ceil
637+
BigDecimal(1).sub(erf(x, high_prec), prec)
638+
end
639+
640+
# Calculates erf(x + a)
641+
private_class_method def _erf_taylor(x, a, erf_a, prec) # :nodoc:
642+
return erf_a if x.zero?
643+
# Let f(x+a) = erf(x+a)*exp((x+a)**2)*sqrt(pi)/2
644+
# = c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ...
645+
# f'(x+a) = 1+2*(x+a)*f(x+a)
646+
# f'(x+a) = c1 + 2*c2*x + 3*c3*x**2 + 4*c4*x**3 + 5*c5*x**4 + ...
647+
# = 1+2*(x+a)*(c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ...)
648+
# therefore,
649+
# c0 = f(a)
650+
# c1 = 2 * a * c0 + 1
651+
# c2 = (2 * c0 + 2 * a * c1) / 2
652+
# c3 = (2 * c1 + 2 * a * c2) / 3
653+
# c4 = (2 * c2 + 2 * a * c3) / 4
654+
#
655+
# All coefficients are positive when a >= 0
656+
657+
scale = BigDecimal(2).div(sqrt(PI(prec), prec), prec)
658+
c_prev = erf_a.div(scale.mult(exp(-a*a, prec), prec), prec)
659+
c_next = (2 * a * c_prev).add(1, prec).mult(x, prec)
660+
sum = c_prev.add(c_next, prec)
661+
662+
2.step do |k|
663+
c = (c_prev.mult(x, prec) + a * c_next).mult(2, prec).mult(x, prec).div(k, prec)
664+
sum = sum.add(c, prec)
665+
c_prev, c_next = c_next, c
666+
break if [c_prev, c_next].all? { |c| c.zero? || (c.exponent < sum.exponent - prec) }
667+
end
668+
value = sum.mult(scale.mult(exp(-(x + a).mult(x + a, prec), prec), prec), prec)
669+
value > 1 ? BigDecimal(1) : value
670+
end
671+
672+
private_class_method def _erfc_asymptotic(x, prec) # :nodoc:
673+
# Let f(x) = erfc(x)*sqrt(pi)*exp(x**2)/2
674+
# f(x) satisfies the following differential equation:
675+
# 2*x*f(x) = f'(x) + 1
676+
# From the above equation, we can derive the following asymptotic expansion:
677+
# f(x) = (0..kmax).sum { (-1)**k * (2*k)! / 4**k / k! / x**(2*k)) } / x
678+
679+
# This asymptotic expansion does not converge.
680+
# But if there is a k that satisfies (2*k)! / 4**k / k! / x**(2*k) < 10**(-prec),
681+
# It is enough to calculate erfc within the given precision.
682+
# Using Stirling's approximation, we can simplify this condition to:
683+
# sqrt(2)/2 + k*log(k) - k - 2*k*log(x) < -prec*log(10)
684+
# and the left side is minimized when k = x**2.
685+
prec += BigDecimal.double_fig
686+
xf = x.to_f
687+
kmax = (1..(xf ** 2).floor).bsearch do |k|
688+
Math.log(2) / 2 + k * Math.log(k) - k - 2 * k * Math.log(xf) < -prec * Math.log(10)
689+
end
690+
return unless kmax
691+
692+
sum = BigDecimal(1)
693+
x2 = x.mult(x, prec)
694+
d = BigDecimal(1)
695+
(1..kmax).each do |k|
696+
d = d.div(x2, prec).mult(1 - 2 * k, prec).div(2, prec)
697+
sum = sum.add(d, prec)
698+
end
699+
sum.div(exp(x2, prec).mult(PI(prec).sqrt(prec), prec), prec).div(x, prec)
700+
end
571701

572702
# call-seq:
573703
# PI(numeric) -> BigDecimal

‎test/bigdecimal/test_bigmath.rb‎

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,8 @@ def test_consistent_precision_acceptance
8282
assert_consistent_precision_acceptance {|prec| BigMath.log10(x, prec) }
8383
assert_consistent_precision_acceptance {|prec| BigMath.log1p(x, prec) }
8484
assert_consistent_precision_acceptance {|prec| BigMath.expm1(x, prec) }
85+
assert_consistent_precision_acceptance {|prec| BigMath.erf(x, prec) }
86+
assert_consistent_precision_acceptance {|prec| BigMath.erfc(x, prec) }
8587

8688
assert_consistent_precision_acceptance {|prec| BigMath.E(prec) }
8789
assert_consistent_precision_acceptance {|prec| BigMath.PI(prec) }
@@ -112,6 +114,8 @@ def test_coerce_argument
112114
assert_equal(log10(bd, N), log10(f, N))
113115
assert_equal(log1p(bd, N), log1p(f, N))
114116
assert_equal(expm1(bd, N), expm1(f, N))
117+
assert_equal(erf(bd, N), erf(f, N))
118+
assert_equal(erfc(bd, N), erfc(f, N))
115119
end
116120

117121
def test_sqrt
@@ -469,4 +473,58 @@ def test_expm1
469473
assert_in_exact_precision(exp(BigDecimal("1.23e-10"), 120) - 1, expm1(BigDecimal("1.23e-10"), 100), 100)
470474
assert_in_exact_precision(exp(123, 120) - 1, expm1(BigDecimal("123"), 100), 100)
471475
end
476+
477+
def test_erf
478+
[-0.5, 0.1, 0.3, 2.1, 3.3].each do |x|
479+
assert_in_epsilon(Math.erf(x), BigMath.erf(BigDecimal(x.to_s), N))
480+
end
481+
assert_equal(1, BigMath.erf(PINF, N))
482+
assert_equal(-1, BigMath.erf(MINF, N))
483+
assert_equal(1, BigMath.erf(BigDecimal(1000), 100))
484+
assert_equal(-1, BigMath.erf(BigDecimal(-1000), 100))
485+
assert_not_equal(1, BigMath.erf(BigDecimal(10), 45))
486+
assert_not_equal(1, BigMath.erf(BigDecimal(15), 100))
487+
assert_equal(1, BigMath.erf(BigDecimal('1e400'), 10))
488+
assert_equal(-1, BigMath.erf(BigDecimal('-1e400'), 10))
489+
assert_equal(BigMath.erf(BigDecimal('1e-300'), N) * BigDecimal('1e-100'), BigMath.erf(BigDecimal('1e-400'), N))
490+
assert_equal(
491+
BigDecimal("0.9953222650189527341620692563672529286108917970400600767383523262004372807199951773676290080196806805"),
492+
BigMath.erf(BigDecimal("2"), 100)
493+
)
494+
assert_converge_in_precision {|n| BigMath.erf(BigDecimal("1e-30"), n) }
495+
assert_converge_in_precision {|n| BigMath.erf(BigDecimal("0.3"), n) }
496+
assert_converge_in_precision {|n| BigMath.erf(SQRT2, n) }
497+
end
498+
499+
def test_erfc
500+
[-0.5, 0.1, 0.3, 2.1, 3.3].each do |x|
501+
assert_in_epsilon(Math.erfc(x), BigMath.erfc(BigDecimal(x.to_s), N))
502+
end
503+
assert_equal(0, BigMath.erfc(PINF, N))
504+
assert_equal(2, BigMath.erfc(MINF, N))
505+
assert_equal(0, BigMath.erfc(BigDecimal('1e400'), 10))
506+
assert_equal(2, BigMath.erfc(BigDecimal('-1e400'), 10))
507+
assert_equal(1, BigMath.erfc(BigDecimal('1e-400'), N))
508+
509+
# erfc with taylor series
510+
assert_equal(
511+
BigDecimal("2.088487583762544757000786294957788611560818119321163727012213713938174695833440290610766384285723554e-45"),
512+
BigMath.erfc(BigDecimal("10"), 100)
513+
)
514+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(0.3), n) }
515+
assert_converge_in_precision {|n| BigMath.erfc(SQRT2, n) }
516+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(8), n) }
517+
# erfc with asymptotic expansion
518+
assert_equal(
519+
BigDecimal("1.896961059966276509268278259713415434936907563929186183462834752900411805205111886605256690776760041e-697"),
520+
BigMath.erfc(BigDecimal("40"), 100)
521+
)
522+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(30), n) }
523+
assert_converge_in_precision {|n| BigMath.erfc(30 * SQRT2, n) }
524+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(50), n) }
525+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(60000), n) }
526+
# Near crossover point between taylor series and asymptotic expansion around prec=150
527+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(19.5), n) }
528+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(20.5), n) }
529+
end
472530
end

‎test/bigdecimal/test_jruby.rb‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,8 @@ def test_bigmath
7171
assert_in_delta(Math.log10(3), BigMath.log10(BigDecimal(3), N))
7272
assert_in_delta(Math.log1p(0.1), BigMath.log1p(BigDecimal('0.1'), N)) if defined? Math.log1p
7373
assert_in_delta(Math.expm1(0.1), BigMath.expm1(BigDecimal('0.1'), N)) if defined? Math.expm1
74+
assert_in_delta(Math.erf(1), BigMath.erf(BigDecimal(1), N))
75+
assert_in_delta(Math.erfc(10), BigMath.erfc(BigDecimal(10), N))
7476
assert_in_delta(Math::PI, BigMath.PI(N))
7577
assert_in_delta(Math::E, BigMath.E(N))
7678
end

0 commit comments

Comments
 (0)