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

Uses asymptotic expansion of erfc if possible and fallback to taylor series of erf

Author: tompng

lib/bigdecimal/math.rb

@@ -8,6 +8,8 @@
 #   sin (x, prec)
 #   cos (x, prec)
 #   atan(x, prec)
+#   erf (x, prec)
+#   erfc(x, prec)
 #   PI  (prec)
 #   E   (prec) == exp(1.0,prec)
 #
@@ -231,4 +233,133 @@ def E(prec)
     raise ArgumentError, "Zero or negative precision for E" if prec <= 0
     BigMath.exp(1, prec)
   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'), 16).to_s
+  #   #=> "0.84270079294971486934122063508259e0"
+  #
+  def erf(x, prec)
+    raise ArgumentError, "Zero or negative precision for erf" if prec <= 0
+    return BigDecimal("NaN") if x.nan?
+    return BigDecimal(0) if x == 0
+    return -erf(-x, prec) if x < 0
+
+    prec += BigDecimal.double_fig
+
+    if x > 8 && (erfc1 = _erfc_asymptotic(x.abs, 1))
+      erfc2 = _erfc_asymptotic(x.abs, [prec + erfc1.exponent, 1].max)
+      return BigDecimal(1).sub(erfc2, prec) if erfc2
+    end
+
+    base = BigDecimal::BASE ** 2
+    x_smallprec = (x * base).fix / base
+    # Taylor series of x with small precision is fast
+    erf1 = _erf_taylor(x_smallprec, BigDecimal(0), BigDecimal(0), prec)
+    # Taylor series converges quickly for small x
+    v = _erf_taylor(x - x_smallprec, x_smallprec, erf1, prec)
+    [BigDecimal(1), v].min
+  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'), 16).to_s
+  #   #=> "0.20884875837625447570007862949578e-44"
+  #
+  def erfc(x, prec)
+    raise ArgumentError, "Zero or negative precision for erfc" if prec <= 0
+    return BigDecimal("NaN") if x.nan?
+    return BigDecimal(1).sub(erf(x, prec), prec + BigDecimal.double_fig) if x < 0
+
+    if x >= 8
+      y = _erfc_asymptotic(x, prec)
+      return y if y
+    end
+
+    prec += BigDecimal.double_fig
+
+    # 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 + (x**2 / log10).ceil
+    BigDecimal(1).sub(erf(x, high_prec), prec)
+  end
+
+
+  private def _erf_taylor(x, a, erf_a, prec)
+    # Let f(x) = erf(x+a)*exp((x+a)**2)*sqrt(pi)/2
+    #          = c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ...
+    # f'(x) = 1+2*(x+a)*f(x)
+    # f'(x) = 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(0)
+    # 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
+
+    return erf_a if x.zero?
+
+    scale = BigDecimal(2).div(sqrt(PI(prec), prec), prec)
+    c_prev = erf_a.div(scale.mult(BigMath.exp(-a*a, prec), prec), prec)
+    c_next = (2 * a * c_prev).add(1, prec).mult(x, prec)
+    v = 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)
+      v = v.add(c, prec)
+      c_prev, c_next = c_next, c
+      break if [c_prev, c_next].all? { |c| c.zero?  || (c.exponent < v.exponent - prec) }
+    end
+    v = v.mult(scale.mult(BigMath.exp(-(x + a).mult(x + a, prec), prec), prec), prec)
+    v > 1 ? BigDecimal(1) : v
+  end
+
+  private def _erfc_asymptotic(x, prec)
+    # 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) = 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.
+    # (2*k)! / 4**k / k! can be approximated as sqrt(2) * (k/e)**k by using Stirling's approximation.
+    prec += BigDecimal.double_fig
+    xf = x.to_f
+    log10xf = Math.log10(xf)
+    kmax = 1
+    until kmax * Math.log10(kmax / Math::E) + 1 - 2 * kmax * log10xf < -prec
+      kmax += 1
+      return if xf * xf < kmax # Unable to calculate with the given precision
+    end
+
+    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
+    expx2 = BigMath.exp(x.mult(x, prec), prec)
+
+    # Workaround for https://github.com/ruby/bigdecimal/issues/345
+    expx2 = BigDecimal(expx2) unless expx2.is_a?(BigDecimal)
+
+    sum.div(expx2.mult(PI(prec).sqrt(prec), prec), prec).div(x, prec)
+  end
 end

test/bigdecimal/test_bigmath.rb

@@ -130,4 +130,45 @@ def test_log
     end
     SRC
   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(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_in_epsilon(
+      BigDecimal("0.995322265018952734162069256367252928610891797040060076738352326200437280719995177367629008019680680487939328715594755785"),
+      BigMath.erf(BigDecimal("2"), 100),
+      BigDecimal("1e-100")
+    )
+    assert_relative_precision {|n| BigMath.erf(BigDecimal("1e-30"), n) }
+    assert_relative_precision {|n| BigMath.erf(BigDecimal("0.3"), n) }
+    assert_relative_precision {|n| BigMath.erf(BigDecimal("1.2"), 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
+    # erfc with taylor series
+    assert_in_epsilon(
+      BigDecimal("2.08848758376254475700078629495778861156081811932116372701221371393817469583344029061076638428572355398152593923652403986e-45"),
+      BigMath.erfc(BigDecimal("10"), 100),
+      BigDecimal("1e-100")
+    )
+    assert_relative_precision {|n| BigMath.erfc(BigDecimal("0.3"), n) }
+    assert_relative_precision {|n| BigMath.erfc(BigDecimal("1.2"), n) }
+    assert_relative_precision {|n| BigMath.erfc(BigDecimal("30"), n) }
+    # erfc with asymptotic expansion
+    assert_in_epsilon(
+      BigDecimal("1.89696105996627650926827825971341543493690756392918618346283475290041180520511188660525669077676004136530598303468056210e-697"),
+      BigMath.erfc(BigDecimal("40"), 100),
+      BigDecimal("1e-100")
+    )
+    assert_relative_precision {|n| BigMath.erfc(BigDecimal("30"), n) }
+    assert_relative_precision {|n| BigMath.erfc(BigDecimal("50"), n) }
+  end
 end