libgo/go/math/erfinv.go

   1 // Copyright 2017 The Go Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package math
   6
   7 /*
   8         Inverse of the floating-point error function.
   9 */
  10
  11 // This implementation is based on the rational approximation
  12 // of percentage points of normal distribution available from
  13 // http://www.jstor.org/stable/2347330.
  14
  15 const (
  16         // Coefficients for approximation to erf in |x| <= 0.85
  17         a0 = 1.1975323115670912564578e0
  18         a1 = 4.7072688112383978012285e1
  19         a2 = 6.9706266534389598238465e2
  20         a3 = 4.8548868893843886794648e3
  21         a4 = 1.6235862515167575384252e4
  22         a5 = 2.3782041382114385731252e4
  23         a6 = 1.1819493347062294404278e4
  24         a7 = 8.8709406962545514830200e2
  25         b0 = 1.0000000000000000000e0
  26         b1 = 4.2313330701600911252e1
  27         b2 = 6.8718700749205790830e2
  28         b3 = 5.3941960214247511077e3
  29         b4 = 2.1213794301586595867e4
  30         b5 = 3.9307895800092710610e4
  31         b6 = 2.8729085735721942674e4
  32         b7 = 5.2264952788528545610e3
  33         // Coefficients for approximation to erf in 0.85 < |x| <= 1-2*exp(-25)
  34         c0 = 1.42343711074968357734e0
  35         c1 = 4.63033784615654529590e0
  36         c2 = 5.76949722146069140550e0
  37         c3 = 3.64784832476320460504e0
  38         c4 = 1.27045825245236838258e0
  39         c5 = 2.41780725177450611770e-1
  40         c6 = 2.27238449892691845833e-2
  41         c7 = 7.74545014278341407640e-4
  42         d0 = 1.4142135623730950488016887e0
  43         d1 = 2.9036514445419946173133295e0
  44         d2 = 2.3707661626024532365971225e0
  45         d3 = 9.7547832001787427186894837e-1
  46         d4 = 2.0945065210512749128288442e-1
  47         d5 = 2.1494160384252876777097297e-2
  48         d6 = 7.7441459065157709165577218e-4
  49         d7 = 1.4859850019840355905497876e-9
  50         // Coefficients for approximation to erf in 1-2*exp(-25) < |x| < 1
  51         e0 = 6.65790464350110377720e0
  52         e1 = 5.46378491116411436990e0
  53         e2 = 1.78482653991729133580e0
  54         e3 = 2.96560571828504891230e-1
  55         e4 = 2.65321895265761230930e-2
  56         e5 = 1.24266094738807843860e-3
  57         e6 = 2.71155556874348757815e-5
  58         e7 = 2.01033439929228813265e-7
  59         f0 = 1.414213562373095048801689e0
  60         f1 = 8.482908416595164588112026e-1
  61         f2 = 1.936480946950659106176712e-1
  62         f3 = 2.103693768272068968719679e-2
  63         f4 = 1.112800997078859844711555e-3
  64         f5 = 2.611088405080593625138020e-5
  65         f6 = 2.010321207683943062279931e-7
  66         f7 = 2.891024605872965461538222e-15
  67 )
  68
  69 // Erfinv returns the inverse error function of x.
  70 //
  71 // Special cases are:
  72 //      Erfinv(1) = +Inf
  73 //      Erfinv(-1) = -Inf
  74 //      Erfinv(x) = NaN if x < -1 or x > 1
  75 //      Erfinv(NaN) = NaN
  76 func Erfinv(x float64) float64 {
  77         // special cases
  78         if IsNaN(x) || x <= -1 || x >= 1 {
  79                 if x == -1 || x == 1 {
  80                         return Inf(int(x))
  81                 }
  82                 return NaN()
  83         }
  84
  85         sign := false
  86         if x < 0 {
  87                 x = -x
  88                 sign = true
  89         }
  90
  91         var ans float64
  92         if x <= 0.85 { // |x| <= 0.85
  93                 r := 0.180625 - 0.25*x*x
  94                 z1 := ((((((a7*r+a6)*r+a5)*r+a4)*r+a3)*r+a2)*r+a1)*r + a0
  95                 z2 := ((((((b7*r+b6)*r+b5)*r+b4)*r+b3)*r+b2)*r+b1)*r + b0
  96                 ans = (x * z1) / z2
  97         } else {
  98                 var z1, z2 float64
  99                 r := Sqrt(Ln2 - Log(1.0-x))
 100                 if r <= 5.0 {
 101                         r -= 1.6
 102                         z1 = ((((((c7*r+c6)*r+c5)*r+c4)*r+c3)*r+c2)*r+c1)*r + c0
 103                         z2 = ((((((d7*r+d6)*r+d5)*r+d4)*r+d3)*r+d2)*r+d1)*r + d0
 104                 } else {
 105                         r -= 5.0
 106                         z1 = ((((((e7*r+e6)*r+e5)*r+e4)*r+e3)*r+e2)*r+e1)*r + e0
 107                         z2 = ((((((f7*r+f6)*r+f5)*r+f4)*r+f3)*r+f2)*r+f1)*r + f0
 108                 }
 109                 ans = z1 / z2
 110         }
 111
 112         if sign {
 113                 return -ans
 114         }
 115         return ans
 116 }
 117
 118 // Erfcinv returns the inverse of Erfc(x).
 119 //
 120 // Special cases are:
 121 //      Erfcinv(0) = +Inf
 122 //      Erfcinv(2) = -Inf
 123 //      Erfcinv(x) = NaN if x < 0 or x > 2
 124 //      Erfcinv(NaN) = NaN
 125 func Erfcinv(x float64) float64 {
 126         return Erfinv(1 - x)
 127 }