gamma.go

Documentation: math

		 1  // Copyright 2010 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  // The original C code, the long comment, and the constants
		 8  // below are from http://netlib.sandia.gov/cephes/cprob/gamma.c.
		 9  // The go code is a simplified version of the original C.
		10  //
		11  //			tgamma.c
		12  //
		13  //			Gamma function
		14  //
		15  // SYNOPSIS:
		16  //
		17  // double x, y, tgamma();
		18  // extern int signgam;
		19  //
		20  // y = tgamma( x );
		21  //
		22  // DESCRIPTION:
		23  //
		24  // Returns gamma function of the argument. The result is
		25  // correctly signed, and the sign (+1 or -1) is also
		26  // returned in a global (extern) variable named signgam.
		27  // This variable is also filled in by the logarithmic gamma
		28  // function lgamma().
		29  //
		30  // Arguments |x| <= 34 are reduced by recurrence and the function
		31  // approximated by a rational function of degree 6/7 in the
		32  // interval (2,3).	Large arguments are handled by Stirling's
		33  // formula. Large negative arguments are made positive using
		34  // a reflection formula.
		35  //
		36  // ACCURACY:
		37  //
		38  //											Relative error:
		39  // arithmetic	 domain		 # trials			peak				 rms
		40  //		DEC			-34, 34			10000			 1.3e-16		 2.5e-17
		41  //		IEEE		-170,-33			20000			 2.3e-15		 3.3e-16
		42  //		IEEE		 -33,	33		 20000			 9.4e-16		 2.2e-16
		43  //		IEEE			33, 171.6	 20000			 2.3e-15		 3.2e-16
		44  //
		45  // Error for arguments outside the test range will be larger
		46  // owing to error amplification by the exponential function.
		47  //
		48  // Cephes Math Library Release 2.8:	June, 2000
		49  // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
		50  //
		51  // The readme file at http://netlib.sandia.gov/cephes/ says:
		52  //		Some software in this archive may be from the book _Methods and
		53  // Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
		54  // International, 1989) or from the Cephes Mathematical Library, a
		55  // commercial product. In either event, it is copyrighted by the author.
		56  // What you see here may be used freely but it comes with no support or
		57  // guarantee.
		58  //
		59  //	 The two known misprints in the book are repaired here in the
		60  // source listings for the gamma function and the incomplete beta
		61  // integral.
		62  //
		63  //	 Stephen L. Moshier
		64  //	 [email protected]
		65  
		66  var _gamP = [...]float64{
		67  	1.60119522476751861407e-04,
		68  	1.19135147006586384913e-03,
		69  	1.04213797561761569935e-02,
		70  	4.76367800457137231464e-02,
		71  	2.07448227648435975150e-01,
		72  	4.94214826801497100753e-01,
		73  	9.99999999999999996796e-01,
		74  }
		75  var _gamQ = [...]float64{
		76  	-2.31581873324120129819e-05,
		77  	5.39605580493303397842e-04,
		78  	-4.45641913851797240494e-03,
		79  	1.18139785222060435552e-02,
		80  	3.58236398605498653373e-02,
		81  	-2.34591795718243348568e-01,
		82  	7.14304917030273074085e-02,
		83  	1.00000000000000000320e+00,
		84  }
		85  var _gamS = [...]float64{
		86  	7.87311395793093628397e-04,
		87  	-2.29549961613378126380e-04,
		88  	-2.68132617805781232825e-03,
		89  	3.47222221605458667310e-03,
		90  	8.33333333333482257126e-02,
		91  }
		92  
		93  // Gamma function computed by Stirling's formula.
		94  // The pair of results must be multiplied together to get the actual answer.
		95  // The multiplication is left to the caller so that, if careful, the caller can avoid
		96  // infinity for 172 <= x <= 180.
		97  // The polynomial is valid for 33 <= x <= 172; larger values are only used
		98  // in reciprocal and produce denormalized floats. The lower precision there
		99  // masks any imprecision in the polynomial.
	 100  func stirling(x float64) (float64, float64) {
	 101  	if x > 200 {
	 102  		return Inf(1), 1
	 103  	}
	 104  	const (
	 105  		SqrtTwoPi	 = 2.506628274631000502417
	 106  		MaxStirling = 143.01608
	 107  	)
	 108  	w := 1 / x
	 109  	w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4])
	 110  	y1 := Exp(x)
	 111  	y2 := 1.0
	 112  	if x > MaxStirling { // avoid Pow() overflow
	 113  		v := Pow(x, 0.5*x-0.25)
	 114  		y1, y2 = v, v/y1
	 115  	} else {
	 116  		y1 = Pow(x, x-0.5) / y1
	 117  	}
	 118  	return y1, SqrtTwoPi * w * y2
	 119  }
	 120  
	 121  // Gamma returns the Gamma function of x.
	 122  //
	 123  // Special cases are:
	 124  //	Gamma(+Inf) = +Inf
	 125  //	Gamma(+0) = +Inf
	 126  //	Gamma(-0) = -Inf
	 127  //	Gamma(x) = NaN for integer x < 0
	 128  //	Gamma(-Inf) = NaN
	 129  //	Gamma(NaN) = NaN
	 130  func Gamma(x float64) float64 {
	 131  	const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
	 132  	// special cases
	 133  	switch {
	 134  	case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
	 135  		return NaN()
	 136  	case IsInf(x, 1):
	 137  		return Inf(1)
	 138  	case x == 0:
	 139  		if Signbit(x) {
	 140  			return Inf(-1)
	 141  		}
	 142  		return Inf(1)
	 143  	}
	 144  	q := Abs(x)
	 145  	p := Floor(q)
	 146  	if q > 33 {
	 147  		if x >= 0 {
	 148  			y1, y2 := stirling(x)
	 149  			return y1 * y2
	 150  		}
	 151  		// Note: x is negative but (checked above) not a negative integer,
	 152  		// so x must be small enough to be in range for conversion to int64.
	 153  		// If |x| were >= 2⁶³ it would have to be an integer.
	 154  		signgam := 1
	 155  		if ip := int64(p); ip&1 == 0 {
	 156  			signgam = -1
	 157  		}
	 158  		z := q - p
	 159  		if z > 0.5 {
	 160  			p = p + 1
	 161  			z = q - p
	 162  		}
	 163  		z = q * Sin(Pi*z)
	 164  		if z == 0 {
	 165  			return Inf(signgam)
	 166  		}
	 167  		sq1, sq2 := stirling(q)
	 168  		absz := Abs(z)
	 169  		d := absz * sq1 * sq2
	 170  		if IsInf(d, 0) {
	 171  			z = Pi / absz / sq1 / sq2
	 172  		} else {
	 173  			z = Pi / d
	 174  		}
	 175  		return float64(signgam) * z
	 176  	}
	 177  
	 178  	// Reduce argument
	 179  	z := 1.0
	 180  	for x >= 3 {
	 181  		x = x - 1
	 182  		z = z * x
	 183  	}
	 184  	for x < 0 {
	 185  		if x > -1e-09 {
	 186  			goto small
	 187  		}
	 188  		z = z / x
	 189  		x = x + 1
	 190  	}
	 191  	for x < 2 {
	 192  		if x < 1e-09 {
	 193  			goto small
	 194  		}
	 195  		z = z / x
	 196  		x = x + 1
	 197  	}
	 198  
	 199  	if x == 2 {
	 200  		return z
	 201  	}
	 202  
	 203  	x = x - 2
	 204  	p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
	 205  	q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
	 206  	return z * p / q
	 207  
	 208  small:
	 209  	if x == 0 {
	 210  		return Inf(1)
	 211  	}
	 212  	return z / ((1 + Euler*x) * x)
	 213  }
	 214  
	 215  func isNegInt(x float64) bool {
	 216  	if x < 0 {
	 217  		_, xf := Modf(x)
	 218  		return xf == 0
	 219  	}
	 220  	return false
	 221  }
	 222
View as plain text