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 5package 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// moshier@na-net.ornl.gov 65 66var _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} 75var _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} 85var _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 polynomial is valid for 33 <= x <= 172. 95func stirling(x float64) float64 { 96 const ( 97 SqrtTwoPi = 2.506628274631000502417 98 MaxStirling = 143.01608 99 ) 100 w := 1 / x 101 w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4]) 102 y := Exp(x) 103 if x > MaxStirling { // avoid Pow() overflow 104 v := Pow(x, 0.5*x-0.25) 105 y = v * (v / y) 106 } else { 107 y = Pow(x, x-0.5) / y 108 } 109 y = SqrtTwoPi * y * w 110 return y 111} 112 113// Gamma returns the Gamma function of x. 114// 115// Special cases are: 116// Gamma(+Inf) = +Inf 117// Gamma(+0) = +Inf 118// Gamma(-0) = -Inf 119// Gamma(x) = NaN for integer x < 0 120// Gamma(-Inf) = NaN 121// Gamma(NaN) = NaN 122func Gamma(x float64) float64 { 123 const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620 124 // special cases 125 switch { 126 case isNegInt(x) || IsInf(x, -1) || IsNaN(x): 127 return NaN() 128 case x == 0: 129 if Signbit(x) { 130 return Inf(-1) 131 } 132 return Inf(1) 133 case x < -170.5674972726612 || x > 171.61447887182298: 134 return Inf(1) 135 } 136 q := Abs(x) 137 p := Floor(q) 138 if q > 33 { 139 if x >= 0 { 140 return stirling(x) 141 } 142 signgam := 1 143 if ip := int(p); ip&1 == 0 { 144 signgam = -1 145 } 146 z := q - p 147 if z > 0.5 { 148 p = p + 1 149 z = q - p 150 } 151 z = q * Sin(Pi*z) 152 if z == 0 { 153 return Inf(signgam) 154 } 155 z = Pi / (Abs(z) * stirling(q)) 156 return float64(signgam) * z 157 } 158 159 // Reduce argument 160 z := 1.0 161 for x >= 3 { 162 x = x - 1 163 z = z * x 164 } 165 for x < 0 { 166 if x > -1e-09 { 167 goto small 168 } 169 z = z / x 170 x = x + 1 171 } 172 for x < 2 { 173 if x < 1e-09 { 174 goto small 175 } 176 z = z / x 177 x = x + 1 178 } 179 180 if x == 2 { 181 return z 182 } 183 184 x = x - 2 185 p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6] 186 q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7] 187 return z * p / q 188 189small: 190 if x == 0 { 191 return Inf(1) 192 } 193 return z / ((1 + Euler*x) * x) 194} 195 196func isNegInt(x float64) bool { 197 if x < 0 { 198 _, xf := Modf(x) 199 return xf == 0 200 } 201 return false 202} 203