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 // moshier (a] na-net.ornl.gov 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
