Source file src/pkg/math/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	//	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