Cody's Algorithm for ERF

Here's a more substantial example of the use of MathP. Erf(x) is the error function - a function related to the cumulative normal distribution. One way to compute the error function is with Cody's algorithm. I've implemented it in Common Lisp to double precision with MathP.

The implementation here starts by defining two constants (i.e. constants in the following equations), then three helper functions. The helper functions are the three rational functions denoted Rlm in Cody's paper. Here's the code in one big MathP expression:

#M{e = exp(1d0) irpi = 1/sqrt(pi); irp = inverse root pi
rlmErf(x) = (3.209377589138469472562d3;; The 4th degree polynomial from Table II.
      + x * (3.774852376853020208137d2
      + x * (1.138641541510501556495d2
      + x * (3.161123743870565596947d0
      + x *  1.857777061846031526730d\-1))))
      /     (2.844236833439170622273d3
      + x * (1.282616526077372275645d3
      + x * (2.440246379344441733056d2
      + x * (2.360129095234412093499d1 + x))))
rlmErc(x) = (1.23033935479799725272d3;; The 8th degree polynomial from Table III.
      + x * (2.05107837782607146532d3
      + x * (1.71204761263407058314d3
      + x * (8.81952221241769090411d2
      + x * (2.98635138197400131132d2
      + x * (6.61191906371416294775d1
      + x * (8.88314979438837594118d0
      + x * (5.64188496988670089180d\-1
      + x *  2.15311535474403846343d\-8))))))))
      /     (1.23033935480374942043d3
      + x * (3.43936767414372163696d3
      + x * (4.36261909014324715820d3
      + x * (3.29079923573345962678d3
      + x * (1.62138957456669018874d3
      + x * (5.37181101862009857509d2
      + x * (1.17693950891312499305d2
      + x * (1.57449261107098347253d1 + x))))))))
rlmEr2(x) = (-6.58749161529837803157d\-4;; The 5th degree polynomial from Table IV.
      + x * (-1.60837851487422766278d\-2
      + x * (-1.25781726111229246204d\-1
      + x * (-3.60344899949804439429d\-1
      + x * (-3.05326634961232344035d\-1
      + x *  -1.63153871373020978498d\-2)))))
      /     (2.33520497626869185443d\-3
      + x * (6.05183413124413191178d\-2
      + x * (5.27905102951428412248d\-1
      + x * (1.87295284992346047209d0
      + x * (2.56852019228982242072d0 + x)))))
defun erfc(x) {"Placeholder" x}
defun erf (x) {"Error function computed with Cody's algorithm"
    x = float(x pi); Make sure it's a double
    if x < 0.0 (-erf(-x))
    if x < 0.5 x * rlmErf(x^2)
               1d0 - erfc(x)}
defun erfc(x) {"Complementary error function computed with Cody's algorithm"
    x = float(x pi); Make sure it's a double
    if x < 0.5 1d0+erf(-x); Based on identities.
    if x < 4.0 e^-x^2 * rlmErc(x);; These two avoid subtraction from one.
               (oox2=1d0/x^2  e^-x^2/x * {irpi + oox2*rlmEr2(oox2)})}}

To check that it works, we have here a small selection of values for erf and erfc, sampled from the Wikipedia page, compared to those calculated using the above code:

erf(0.1) ≈ 0.1124629 ≈ 0.112463d0
erfc(15)/2 ≈ 3.6065e-100 ≈ 3.606497d-100
erf(0.95) ≈ 0.8208908 ≈ 0.820891d0
erfc(0.7) ≈ 0.3221988 ≈ 0.322199d0

Going a bit further, here are some values based on the arbitrarily precise values computed by Wolfram Alpha:

erf(1/100)
0.0112834155558496160
0.0112834155558496169...
erf(1/3)
0.3626481117660628000
0.3626481117660629334...
erf(1/2)
0.5204998778130465000
0.52049987781304653768...
erf(1)
0.8427007929497148000
0.8427007929497148693...
erf(4)
0.9999999845827421000
0.9999999845827420997...

So it appears (from this small sample at least) that this implementation is accurate to about 1e-16 (see the value for erf(1/3)), but usually 1e-17, which appears to be the resolution of doubles for Corman Common Lisp.