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	//   [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 polynomial is valid for 33 <= x <= 172.
    95	func 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(x) returns the Gamma function of x.
   114	//
   115	// Special cases are:
   116	//	Gamma(±Inf) = ±Inf
   117	//	Gamma(NaN) = NaN
   118	// Large values overflow to +Inf.
   119	// Zero and negative integer arguments return ±Inf.
   120	func Gamma(x float64) float64 {
   121		const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
   122		// special cases
   123		switch {
   124		case IsInf(x, -1) || IsNaN(x):
   125			return x
   126		case x < -170.5674972726612 || x > 171.61447887182298:
   127			return Inf(1)
   128		}
   129		q := Abs(x)
   130		p := Floor(q)
   131		if q > 33 {
   132			if x >= 0 {
   133				return stirling(x)
   134			}
   135			signgam := 1
   136			if ip := int(p); ip&1 == 0 {
   137				signgam = -1
   138			}
   139			z := q - p
   140			if z > 0.5 {
   141				p = p + 1
   142				z = q - p
   143			}
   144			z = q * Sin(Pi*z)
   145			if z == 0 {
   146				return Inf(signgam)
   147			}
   148			z = Pi / (Abs(z) * stirling(q))
   149			return float64(signgam) * z
   150		}
   151	
   152		// Reduce argument
   153		z := 1.0
   154		for x >= 3 {
   155			x = x - 1
   156			z = z * x
   157		}
   158		for x < 0 {
   159			if x > -1e-09 {
   160				goto small
   161			}
   162			z = z / x
   163			x = x + 1
   164		}
   165		for x < 2 {
   166			if x < 1e-09 {
   167				goto small
   168			}
   169			z = z / x
   170			x = x + 1
   171		}
   172	
   173		if x == 2 {
   174			return z
   175		}
   176	
   177		x = x - 2
   178		p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
   179		q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
   180		return z * p / q
   181	
   182	small:
   183		if x == 0 {
   184			return Inf(1)
   185		}
   186		return z / ((1 + Euler*x) * x)
   187	}