Source file src/pkg/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	// https://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	}
   128	

View as plain text