Skip to content

Commit 75582a4

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 84646b8 commit 75582a4

File tree

2 files changed

+180
-0
lines changed

2 files changed

+180
-0
lines changed

‎lib/bigdecimal/math.rb‎

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

572700
# call-seq:
573701
# PI(numeric) -> BigDecimal

‎test/bigdecimal/test_bigmath.rb‎

Lines changed: 52 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,52 @@ 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(
488+
BigDecimal("0.9953222650189527341620692563672529286108917970400600767383523262004372807199951773676290080196806805"),
489+
BigMath.erf(BigDecimal("2"), 100)
490+
)
491+
assert_converge_in_precision {|n| BigMath.erf(BigDecimal("1e-30"), n) }
492+
assert_converge_in_precision {|n| BigMath.erf(BigDecimal("0.3"), n) }
493+
assert_converge_in_precision {|n| BigMath.erf(SQRT2, n) }
494+
end
495+
496+
def test_erfc
497+
[-0.5, 0.1, 0.3, 2.1, 3.3].each do |x|
498+
assert_in_epsilon(Math.erfc(x), BigMath.erfc(BigDecimal(x.to_s), N))
499+
end
500+
assert_equal(0, BigMath.erfc(PINF, N))
501+
assert_equal(2, BigMath.erfc(MINF, N))
502+
503+
# erfc with taylor series
504+
assert_equal(
505+
BigDecimal("2.088487583762544757000786294957788611560818119321163727012213713938174695833440290610766384285723554e-45"),
506+
BigMath.erfc(BigDecimal("10"), 100)
507+
)
508+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(0.3), n) }
509+
assert_converge_in_precision {|n| BigMath.erfc(SQRT2, n) }
510+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(8), n) }
511+
# erfc with asymptotic expansion
512+
assert_equal(
513+
BigDecimal("1.896961059966276509268278259713415434936907563929186183462834752900411805205111886605256690776760041e-697"),
514+
BigMath.erfc(BigDecimal("40"), 100)
515+
)
516+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(30), n) }
517+
assert_converge_in_precision {|n| BigMath.erfc(30 * SQRT2, n) }
518+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(50), n) }
519+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(60000), n) }
520+
# Near crossover point between taylor series and asymptotic expansion around prec=150
521+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(19.5), n) }
522+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(20.5), n) }
523+
end
472524
end

0 commit comments

Comments
 (0)