Skip to content

Commit df97ecb

Browse files
mrknmrkn
authored
Merge branch 'master' into bigmath_gamma
2 parents 174e1a0 + d1a2bc1 commit df97ecb

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
# gamma(x, prec)
2830
# lgamma(x, prec)
2931
# PI (prec)
@@ -570,6 +572,134 @@ def expm1(x, prec)
570572
exp_prec > 0 ? exp(x, exp_prec).sub(1, prec) : BigDecimal(-1)
571573
end
572574

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

‎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
assert_consistent_precision_acceptance {|prec| BigMath.gamma(x, prec) }
8688
assert_consistent_precision_acceptance {|prec| BigMath.lgamma(x, prec) }
8789

@@ -114,6 +116,8 @@ def test_coerce_argument
114116
assert_equal(log10(bd, N), log10(f, N))
115117
assert_equal(log1p(bd, N), log1p(f, N))
116118
assert_equal(expm1(bd, N), expm1(f, N))
119+
assert_equal(erf(bd, N), erf(f, N))
120+
assert_equal(erfc(bd, N), erfc(f, N))
117121
assert_equal(gamma(bd, N), gamma(f, N))
118122
assert_equal(lgamma(bd, N), lgamma(f, N))
119123
end
@@ -474,6 +478,60 @@ def test_expm1
474478
assert_in_exact_precision(exp(123, 120) - 1, expm1(BigDecimal("123"), 100), 100)
475479
end
476480

481+
def test_erf
482+
[-0.5, 0.1, 0.3, 2.1, 3.3].each do |x|
483+
assert_in_epsilon(Math.erf(x), BigMath.erf(BigDecimal(x.to_s), N))
484+
end
485+
assert_equal(1, BigMath.erf(PINF, N))
486+
assert_equal(-1, BigMath.erf(MINF, N))
487+
assert_equal(1, BigMath.erf(BigDecimal(1000), 100))
488+
assert_equal(-1, BigMath.erf(BigDecimal(-1000), 100))
489+
assert_not_equal(1, BigMath.erf(BigDecimal(10), 45))
490+
assert_not_equal(1, BigMath.erf(BigDecimal(15), 100))
491+
assert_equal(1, BigMath.erf(BigDecimal('1e400'), 10))
492+
assert_equal(-1, BigMath.erf(BigDecimal('-1e400'), 10))
493+
assert_equal(BigMath.erf(BigDecimal('1e-300'), N) * BigDecimal('1e-100'), BigMath.erf(BigDecimal('1e-400'), N))
494+
assert_equal(
495+
BigDecimal("0.9953222650189527341620692563672529286108917970400600767383523262004372807199951773676290080196806805"),
496+
BigMath.erf(BigDecimal("2"), 100)
497+
)
498+
assert_converge_in_precision {|n| BigMath.erf(BigDecimal("1e-30"), n) }
499+
assert_converge_in_precision {|n| BigMath.erf(BigDecimal("0.3"), n) }
500+
assert_converge_in_precision {|n| BigMath.erf(SQRT2, n) }
501+
end
502+
503+
def test_erfc
504+
[-0.5, 0.1, 0.3, 2.1, 3.3].each do |x|
505+
assert_in_epsilon(Math.erfc(x), BigMath.erfc(BigDecimal(x.to_s), N))
506+
end
507+
assert_equal(0, BigMath.erfc(PINF, N))
508+
assert_equal(2, BigMath.erfc(MINF, N))
509+
assert_equal(0, BigMath.erfc(BigDecimal('1e400'), 10))
510+
assert_equal(2, BigMath.erfc(BigDecimal('-1e400'), 10))
511+
assert_equal(1, BigMath.erfc(BigDecimal('1e-400'), N))
512+
513+
# erfc with taylor series
514+
assert_equal(
515+
BigDecimal("2.088487583762544757000786294957788611560818119321163727012213713938174695833440290610766384285723554e-45"),
516+
BigMath.erfc(BigDecimal("10"), 100)
517+
)
518+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(0.3), n) }
519+
assert_converge_in_precision {|n| BigMath.erfc(SQRT2, n) }
520+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(8), n) }
521+
# erfc with asymptotic expansion
522+
assert_equal(
523+
BigDecimal("1.896961059966276509268278259713415434936907563929186183462834752900411805205111886605256690776760041e-697"),
524+
BigMath.erfc(BigDecimal("40"), 100)
525+
)
526+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(30), n) }
527+
assert_converge_in_precision {|n| BigMath.erfc(30 * SQRT2, n) }
528+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(50), n) }
529+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(60000), n) }
530+
# Near crossover point between taylor series and asymptotic expansion around prec=150
531+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(19.5), n) }
532+
assert_converge_in_precision {|n| BigMath.erfc(BigDecimal(20.5), n) }
533+
end
534+
477535
def test_gamma
478536
[-1.8, -0.7, 0.6, 1.5, 2.4].each do |x|
479537
assert_in_epsilon(Math.gamma(x), gamma(BigDecimal(x.to_s), N))

‎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.gamma(0.5), BigMath.gamma(BigDecimal('0.5'), N))
7577
assert_in_delta(Math.lgamma(0.5).first, BigMath.lgamma(BigDecimal('0.5'), N).first)
7678
assert_in_delta(Math::PI, BigMath.PI(N))

0 commit comments

Comments
 (0)