erfinv.go

Documentation: math

		 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