Merge pull request #2987 from stevengj/erfinv

JeffBezanson · JeffBezanson · commit ccceddc3d80b · 2013-05-01T12:13:21.000-07:00
add inverse error functions erfinv and erfcinv
diff --git a/base/exports.jl b/base/exports.jl
@@ -318,6 +318,8 @@ export
     erfc,
     erfcx,
     erfi,
+    erfinv,
+    erfcinv,
     exp,
     exp2,
     expm1,
diff --git a/base/math.jl b/base/math.jl
@@ -14,7 +14,8 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
        airy, airyai, airyprime, airyaiprime, airybi, airybiprime,
        besselj0, besselj1, besselj, bessely0, bessely1, bessely,
        hankelh1, hankelh2, besseli, besselk, besselh,
-       beta, lbeta, eta, zeta, digamma
+       beta, lbeta, eta, zeta, digamma,
+       erfinv, erfcinv
 
 import Base.log, Base.exp, Base.sin, Base.cos, Base.tan, Base.sinh, Base.cosh,
        Base.tanh, Base.asin, Base.acos, Base.atan, Base.asinh, Base.acosh,
@@ -695,4 +696,214 @@ for f in (:erfcx, :erfi, :Dawson)
     end
 end
 
+# evaluate p[1] + x * (p[2] + x * (....)), i.e. a polynomial via Horner's rule
+macro horner(x, p...)
+    ex = p[end]
+    for i = length(p)-1:-1:1
+        ex = :($(p[i]) + $x * $ex)
+    end
+    ex
+end
+
+# Compute the inverse of the error function: erf(erfinv(x)) == x, 
+# using the rational approximants tabulated in:
+#     J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev 
+#     approximations for the inverse of the error function," Math. Comp. 30,
+#     pp. 827--830 (1976).
+#         http://dx.doi.org/10.1090/S0025-5718-1976-0421040-7 
+#         http://www.jstor.org/stable/2005402
+function erfinv(x::Float64)
+    a = abs(x)
+    if a >= 1.0
+        if x == 1.0
+            return inf(Float64)
+        elseif x == -1.0
+            return -inf(Float64)
+        end
+        throw(DomainError())
+    elseif a <= 0.75 # Table 17 in Blair et al.
+        t = x*x - 0.5625
+        return x * @horner(t, 0.16030_49558_44066_229311e2, 
+                             -0.90784_95926_29603_26650e2,
+                              0.18644_91486_16209_87391e3,
+                             -0.16900_14273_46423_82420e3,
+                              0.65454_66284_79448_7048e2,
+                             -0.86421_30115_87247_794e1,
+                              0.17605_87821_39059_0) /
+                   @horner(t, 0.14780_64707_15138_316110e2,
+                             -0.91374_16702_42603_13936e2,
+                              0.21015_79048_62053_17714e3,
+                             -0.22210_25412_18551_32366e3,
+                              0.10760_45391_60551_23830e3,
+                             -0.20601_07303_28265_443e2,
+                              0.1e1)
+    elseif a <= 0.9375 # Table 37 in Blair et al.
+        t = x*x - 0.87890625
+        return x * @horner(t, -0.15238_92634_40726_128e-1,
+                               0.34445_56924_13612_5216,
+                              -0.29344_39867_25424_78687e1,
+                               0.11763_50570_52178_27302e2,
+                              -0.22655_29282_31011_04193e2,
+                               0.19121_33439_65803_30163e2,
+                              -0.54789_27619_59831_8769e1,
+                               0.23751_66890_24448) /
+                   @horner(t, -0.10846_51696_02059_954e-1,
+                               0.26106_28885_84307_8511,
+                              -0.24068_31810_43937_57995e1,
+                               0.10695_12997_33870_14469e2,
+                              -0.23716_71552_15965_81025e2,
+                               0.24640_15894_39172_84883e2,
+                              -0.10014_37634_97830_70835e2,
+                               0.1e1)
+    else # Table 57 in Blair et al.
+        t = 1.0 / sqrt(-log(1.0 - a))
+        return @horner(t, 0.10501_31152_37334_38116e-3,
+                          0.10532_61131_42333_38164_25e-1,
+                          0.26987_80273_62432_83544_516,
+                          0.23268_69578_89196_90806_414e1,
+                          0.71678_54794_91079_96810_001e1,
+                          0.85475_61182_21678_27825_185e1,
+                          0.68738_08807_35438_39802_913e1,
+                          0.36270_02483_09587_08930_02e1,
+                          0.88606_27392_96515_46814_9) / 
+              (copysign(t, x) *
+               @horner(t, 0.10501_26668_70303_37690e-3,
+                          0.10532_86230_09333_27531_11e-1,
+                          0.27019_86237_37515_54845_553,
+                          0.23501_43639_79702_53259_123e1,
+                          0.76078_02878_58012_77064_351e1,
+                          0.11181_58610_40569_07827_3451e2,
+                          0.11948_78791_84353_96667_8438e2,
+                          0.81922_40974_72699_07893_913e1,
+                          0.40993_87907_63680_15361_45e1,
+                          0.1e1))
+    end
+end
+
+function erfinv(x::Float32)
+    a = abs(x)
+    if a >= 1.0f0
+        if x == 1.0f0
+            return inf(Float64)
+        elseif x == -1.0f0
+            return -inf(Float64)
+        end
+        throw(DomainError())
+    elseif a <= 0.75f0 # Table 10 in Blair et al.
+        t = x*x - 0.5625f0
+        return x * @horner(t, -0.13095_99674_22f2,
+                               0.26785_22576_0f2,
+                              -0.92890_57365f1) /
+                   @horner(t, -0.12074_94262_97f2,
+                               0.30960_61452_9f2,
+                              -0.17149_97799_1f2,
+                               0.1f1)
+    elseif a <= 0.9375f0 # Table 29 in Blair et al.
+        t = x*x - 0.87890625f0
+        return x * @horner(t, -0.12402_56522_1f0,
+                               0.10688_05957_4f1,
+                              -0.19594_55607_8f1,
+                               0.42305_81357f0) /
+                   @horner(t, -0.88276_97997f-1,
+                               0.89007_43359f0,
+                              -0.21757_03119_6f1,
+                               0.1f1)
+    else # Table 50 in Blair et al.
+        t = 1.0f0 / sqrt(-log(1.0f0 - a))
+        return @horner(t, 0.15504_70003_116f0,
+                          0.13827_19649_631f1,
+                          0.69096_93488_87f0,
+                         -0.11280_81391_617f1,
+                          0.68054_42468_25f0,
+                         -0.16444_15679_1f0) /
+              (copysign(t, x) *
+               @horner(t, 0.15502_48498_22f0,
+                          0.13852_28141_995f1,
+                          0.1f1))
+    end
+end
+
+erfinv(x::Integer) = erfinv(float(x))
+@vectorize_1arg Real erfinv
+
+# Inverse complementary error function: use Blair tables for y = 1-x, 
+# exploiting the greater accuracy of y (vs. x) when y is small.
+function erfcinv(y::Float64)
+    if y > 0.0625
+        return erfinv(1.0 - y)
+    elseif y <= 0.0
+        if y == 0.0
+            return inf(Float64)
+        end
+        throw(DomainError())
+    elseif y >= 1e-100 # Table 57 in Blair et al.
+        t = 1.0 / sqrt(-log(y))
+        return @horner(t, 0.10501_31152_37334_38116e-3,
+                          0.10532_61131_42333_38164_25e-1,
+                          0.26987_80273_62432_83544_516,
+                          0.23268_69578_89196_90806_414e1,
+                          0.71678_54794_91079_96810_001e1,
+                          0.85475_61182_21678_27825_185e1,
+                          0.68738_08807_35438_39802_913e1,
+                          0.36270_02483_09587_08930_02e1,
+                          0.88606_27392_96515_46814_9) / 
+              (t *
+               @horner(t, 0.10501_26668_70303_37690e-3,
+                          0.10532_86230_09333_27531_11e-1,
+                          0.27019_86237_37515_54845_553,
+                          0.23501_43639_79702_53259_123e1,
+                          0.76078_02878_58012_77064_351e1,
+                          0.11181_58610_40569_07827_3451e2,
+                          0.11948_78791_84353_96667_8438e2,
+                          0.81922_40974_72699_07893_913e1,
+                          0.40993_87907_63680_15361_45e1,
+                          0.1e1))
+    else # Table 80 in Blair et al.
+        t = 1.0 / sqrt(-log(y))
+        return @horner(t, 0.34654_29858_80863_50177e-9,
+                          0.25084_67920_24075_70520_55e-6,
+                          0.47378_13196_37286_02986_534e-4,
+                          0.31312_60375_97786_96408_3388e-2,
+                          0.77948_76454_41435_36994_854e-1,
+                          0.70045_68123_35816_43868_271e0,
+                          0.18710_42034_21679_31668_683e1,
+                          0.71452_54774_31351_45428_3e0) /
+          (t * @horner(t, 0.34654_29567_31595_11156e-9,
+                          0.25084_69079_75880_27114_87e-6,
+                          0.47379_53129_59749_13536_339e-4,
+                          0.31320_63536_46177_68848_0813e-2,
+                          0.78073_48906_27648_97214_733e-1,
+                          0.70715_04479_95337_58619_993e0,
+                          0.19998_51543_49112_15105_214e1,
+                          0.15072_90269_27316_80008_56e1,
+                          0.1e1))
+    end
+end
+
+function erfcinv(y::Float32)
+    if y > 0.0625f0
+        return erfinv(1.0f0 - y)
+    elseif y <= 0.0f0
+        if y == 0.0f0
+            return inf(Float32)
+        end
+        throw(DomainError())
+    else # Table 50 in Blair et al.
+        t = 1.0f0 / sqrt(-log(y))
+        return @horner(t, 0.15504_70003_116f0,
+                       0.13827_19649_631f1,
+                       0.69096_93488_87f0,
+                       -0.11280_81391_617f1,
+                       0.68054_42468_25f0,
+                       -0.16444_15679_1f0) /
+        (t *
+         @horner(t, 0.15502_48498_22f0,
+                 0.13852_28141_995f1,
+                 0.1f1))
+    end
+end
+
+erfcinv(x::Integer) = erfcinv(float(x))
+@vectorize_1arg Real erfcinv
+
 end # module
diff --git a/doc/helpdb.jl b/doc/helpdb.jl
@@ -2695,6 +2695,22 @@ FloatingPoint
 
 "),
 
+("Mathematical Functions","Base","erfinv","erfinv(x)
+
+
+   Compute the inverse error function of a real \"x\", 
+   so that \\operatorname{erf}(\\operatorname{erfinv}(x)) = x.
+
+"),
+
+("Mathematical Functions","Base","erfcinv","erfcinv(x)
+
+
+   Compute the inverse complementary error function of a real \"x\", 
+   so that \\operatorname{erfc}(\\operatorname{erfcinv}(x)) = x.
+
+"),
+
 ("Mathematical Functions","Base","dawson","dawson(x)
 
 
diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
@@ -1527,6 +1527,16 @@ Mathematical Functions
    Compute the Dawson function (scaled imaginary error function) of ``x``,
    defined by :math:`\frac{\sqrt{\pi}}{2} e^{-x^2} \operatorname{erfi}(x)`.
 
+.. function:: erfinv(x)
+
+   Compute the inverse error function of a real ``x``,
+   defined by :math:`\operatorname{erf}(\operatorname{erfinv}(x)) = x`.
+
+.. function:: erfcinv(x)
+
+   Compute the inverse error complementary function of a real ``x``,
+   defined by :math:`\operatorname{erfc}(\operatorname{erfcinv}(x)) = x`.
+
 .. function:: real(z)
 
    Return the real part of the complex number ``z``
diff --git a/test/math.jl b/test/math.jl
@@ -9,6 +9,8 @@
 @test_approx_eq erfc(1) 0.15729920705028513066
 @test_approx_eq erfcx(1) 0.42758357615580700442
 @test_approx_eq erfi(1) 1.6504257587975428760
+@test_approx_eq erfinv(0.84270079294971486934) 1
+@test_approx_eq erfcinv(0.15729920705028513066) 1
 @test_approx_eq dawson(1) 0.53807950691276841914
 # TODO: complex versions only supported on 64-bit for now
 if WORD_SIZE==64
@@ -18,6 +20,17 @@ if WORD_SIZE==64
     @test_approx_eq erfi(1+2im) -0.011259006028815025076+1.0036063427256517509im
     @test_approx_eq dawson(1+2im) -13.388927316482919244-11.828715103889593303im
 end
+for x in logspace(-200, -0.01)
+    @test_approx_eq_eps erf(erfinv(x)) x 1e-12*x
+    @test_approx_eq_eps erf(erfinv(-x)) -x 1e-12*x
+    @test_approx_eq_eps erfc(erfcinv(2*x)) 2*x 1e-12*x
+    if x > 1e-20
+        xf = float32(x)
+        @test_approx_eq_eps erf(erfinv(xf)) xf 1e-5*xf
+        @test_approx_eq_eps erf(erfinv(-xf)) -xf 1e-5*xf
+        @test_approx_eq_eps erfc(erfcinv(2xf)) 2xf 1e-5*xf
+    end
+end
 
 # airy
 @test_approx_eq airy(1.8) 0.0470362168668458052247

-Original file line number
+Diff line change
     erfc,
     erfcx,
     erfi,
 +    erfinv,
 +    erfcinv,
     exp,
     exp2,
     expm1,