// Copyright 2010 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package math

// The original C code, the long comment, and the constants
// below are from http://netlib.sandia.gov/cephes/cprob/gamma.c.
// The go code is a simplified version of the original C.
//
//      tgamma.c
//
//      Gamma function
//
// SYNOPSIS:
//
// double x, y, tgamma();
// extern int signgam;
//
// y = tgamma( x );
//
// DESCRIPTION:
//
// Returns gamma function of the argument. The result is
// correctly signed, and the sign (+1 or -1) is also
// returned in a global (extern) variable named signgam.
// This variable is also filled in by the logarithmic gamma
// function lgamma().
//
// Arguments |x| <= 34 are reduced by recurrence and the function
// approximated by a rational function of degree 6/7 in the
// interval (2,3).  Large arguments are handled by Stirling's
// formula. Large negative arguments are made positive using
// a reflection formula.
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain     # trials      peak         rms
//    DEC      -34, 34      10000       1.3e-16     2.5e-17
//    IEEE    -170,-33      20000       2.3e-15     3.3e-16
//    IEEE     -33,  33     20000       9.4e-16     2.2e-16
//    IEEE      33, 171.6   20000       2.3e-15     3.2e-16
//
// Error for arguments outside the test range will be larger
// owing to error amplification by the exponential function.
//
// Cephes Math Library Release 2.8:  June, 2000
// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
//
// The readme file at http://netlib.sandia.gov/cephes/ says:
//    Some software in this archive may be from the book _Methods and
// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
// International, 1989) or from the Cephes Mathematical Library, a
// commercial product. In either event, it is copyrighted by the author.
// What you see here may be used freely but it comes with no support or
// guarantee.
//
//   The two known misprints in the book are repaired here in the
// source listings for the gamma function and the incomplete beta
// integral.
//
//   Stephen L. Moshier
//   moshier@na-net.ornl.gov

var _gamP = [...]float64{
	1.60119522476751861407e-04,
	1.19135147006586384913e-03,
	1.04213797561761569935e-02,
	4.76367800457137231464e-02,
	2.07448227648435975150e-01,
	4.94214826801497100753e-01,
	9.99999999999999996796e-01,
}
var _gamQ = [...]float64{
	-2.31581873324120129819e-05,
	5.39605580493303397842e-04,
	-4.45641913851797240494e-03,
	1.18139785222060435552e-02,
	3.58236398605498653373e-02,
	-2.34591795718243348568e-01,
	7.14304917030273074085e-02,
	1.00000000000000000320e+00,
}
var _gamS = [...]float64{
	7.87311395793093628397e-04,
	-2.29549961613378126380e-04,
	-2.68132617805781232825e-03,
	3.47222221605458667310e-03,
	8.33333333333482257126e-02,
}

// Gamma function computed by Stirling's formula.
// The pair of results must be multiplied together to get the actual answer.
// The multiplication is left to the caller so that, if careful, the caller can avoid
// infinity for 172 <= x <= 180.
// The polynomial is valid for 33 <= x <= 172; larger values are only used
// in reciprocal and produce denormalized floats. The lower precision there
// masks any imprecision in the polynomial.
func ( float64) (float64, float64) {
	if  > 200 {
		return Inf(1), 1
	}
	const (
		   = 2.506628274631000502417
		 = 143.01608
	)
	 := 1 / 
	 = 1 + *((((_gamS[0]*+_gamS[1])*+_gamS[2])*+_gamS[3])*+_gamS[4])
	 := Exp()
	 := 1.0
	if  >  { // avoid Pow() overflow
		 := Pow(, 0.5*-0.25)
		,  = , /
	} else {
		 = Pow(, -0.5) / 
	}
	return ,  *  * 
}

// Gamma returns the Gamma function of x.
//
// Special cases are:
//
//	Gamma(+Inf) = +Inf
//	Gamma(+0) = +Inf
//	Gamma(-0) = -Inf
//	Gamma(x) = NaN for integer x < 0
//	Gamma(-Inf) = NaN
//	Gamma(NaN) = NaN
func ( float64) float64 {
	const  = 0.57721566490153286060651209008240243104215933593992 // A001620
	// special cases
	switch {
	case isNegInt() || IsInf(, -1) || IsNaN():
		return NaN()
	case IsInf(, 1):
		return Inf(1)
	case  == 0:
		if Signbit() {
			return Inf(-1)
		}
		return Inf(1)
	}
	 := Abs()
	 := Floor()
	if  > 33 {
		if  >= 0 {
			,  := stirling()
			return  * 
		}
		// Note: x is negative but (checked above) not a negative integer,
		// so x must be small enough to be in range for conversion to int64.
		// If |x| were >= 2⁶³ it would have to be an integer.
		 := 1
		if  := int64(); &1 == 0 {
			 = -1
		}
		 :=  - 
		if  > 0.5 {
			 =  + 1
			 =  - 
		}
		 =  * Sin(Pi*)
		if  == 0 {
			return Inf()
		}
		,  := stirling()
		 := Abs()
		 :=  *  * 
		if IsInf(, 0) {
			 = Pi /  /  / 
		} else {
			 = Pi / 
		}
		return float64() * 
	}

	// Reduce argument
	 := 1.0
	for  >= 3 {
		 =  - 1
		 =  * 
	}
	for  < 0 {
		if  > -1e-09 {
			goto 
		}
		 =  / 
		 =  + 1
	}
	for  < 2 {
		if  < 1e-09 {
			goto 
		}
		 =  / 
		 =  + 1
	}

	if  == 2 {
		return 
	}

	 =  - 2
	 = (((((*_gamP[0]+_gamP[1])*+_gamP[2])*+_gamP[3])*+_gamP[4])*+_gamP[5])* + _gamP[6]
	 = ((((((*_gamQ[0]+_gamQ[1])*+_gamQ[2])*+_gamQ[3])*+_gamQ[4])*+_gamQ[5])*+_gamQ[6])* + _gamQ[7]
	return  *  / 

:
	if  == 0 {
		return Inf(1)
	}
	return  / ((1 + *) * )
}

func ( float64) bool {
	if  < 0 {
		,  := Modf()
		return  == 0
	}
	return false
}