gamma.go

     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@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  //
   125  //	Gamma(+Inf) = +Inf
   126  //	Gamma(+0) = +Inf
   127  //	Gamma(-0) = -Inf
   128  //	Gamma(x) = NaN for integer x < 0
   129  //	Gamma(-Inf) = NaN
   130  //	Gamma(NaN) = NaN
   131  func Gamma(x float64) float64 {
   132  	const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
   133  	// special cases
   134  	switch {
   135  	case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
   136  		return NaN()
   137  	case IsInf(x, 1):
   138  		return Inf(1)
   139  	case x == 0:
   140  		if Signbit(x) {
   141  			return Inf(-1)
   142  		}
   143  		return Inf(1)
   144  	}
   145  	q := Abs(x)
   146  	p := Floor(q)
   147  	if q > 33 {
   148  		if x >= 0 {
   149  			y1, y2 := stirling(x)
   150  			return y1 * y2
   151  		}
   152  		// Note: x is negative but (checked above) not a negative integer,
   153  		// so x must be small enough to be in range for conversion to int64.
   154  		// If |x| were >= 2⁶³ it would have to be an integer.
   155  		signgam := 1
   156  		if ip := int64(p); ip&1 == 0 {
   157  			signgam = -1
   158  		}
   159  		z := q - p
   160  		if z > 0.5 {
   161  			p = p + 1
   162  			z = q - p
   163  		}
   164  		z = q * Sin(Pi*z)
   165  		if z == 0 {
   166  			return Inf(signgam)
   167  		}
   168  		sq1, sq2 := stirling(q)
   169  		absz := Abs(z)
   170  		d := absz * sq1 * sq2
   171  		if IsInf(d, 0) {
   172  			z = Pi / absz / sq1 / sq2
   173  		} else {
   174  			z = Pi / d
   175  		}
   176  		return float64(signgam) * z
   177  	}
   178  
   179  	// Reduce argument
   180  	z := 1.0
   181  	for x >= 3 {
   182  		x = x - 1
   183  		z = z * x
   184  	}
   185  	for x < 0 {
   186  		if x > -1e-09 {
   187  			goto small
   188  		}
   189  		z = z / x
   190  		x = x + 1
   191  	}
   192  	for x < 2 {
   193  		if x < 1e-09 {
   194  			goto small
   195  		}
   196  		z = z / x
   197  		x = x + 1
   198  	}
   199  
   200  	if x == 2 {
   201  		return z
   202  	}
   203  
   204  	x = x - 2
   205  	p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
   206  	q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
   207  	return z * p / q
   208  
   209  small:
   210  	if x == 0 {
   211  		return Inf(1)
   212  	}
   213  	return z / ((1 + Euler*x) * x)
   214  }
   215  
   216  func isNegInt(x float64) bool {
   217  	if x < 0 {
   218  		_, xf := Modf(x)
   219  		return xf == 0
   220  	}
   221  	return false
   222  }
   223
View as plain text