src/pkg/math/lgamma.go - go - Git at Google

 // 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

 /*
 	Floating-point logarithm of the Gamma function.
 */

 // The original C code and the long comment below are
 // from FreeBSD's /usr/src/lib/msun/src/e_lgamma_r.c and
 // came with this notice.  The go code is a simplified
 // version of the original C.
 //
 // ====================================================
 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 //
 // Developed at SunPro, a Sun Microsystems, Inc. business.
 // Permission to use, copy, modify, and distribute this
 // software is freely granted, provided that this notice
 // is preserved.
 // ====================================================
 //
 // __ieee754_lgamma_r(x, signgamp)
 // Reentrant version of the logarithm of the Gamma function
 // with user provided pointer for the sign of Gamma(x).
 //
 // Method:
 //   1. Argument Reduction for 0 < x <= 8
 //      Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
 //      reduce x to a number in [1.5,2.5] by
 //              lgamma(1+s) = log(s) + lgamma(s)
 //      for example,
 //              lgamma(7.3) = log(6.3) + lgamma(6.3)
 //                          = log(6.3*5.3) + lgamma(5.3)
 //                          = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
 //   2. Polynomial approximation of lgamma around its
 //      minimum (ymin=1.461632144968362245) to maintain monotonicity.
 //      On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
 //              Let z = x-ymin;
 //              lgamma(x) = -1.214862905358496078218 + z**2*poly(z)
 //              poly(z) is a 14 degree polynomial.
 //   2. Rational approximation in the primary interval [2,3]
 //      We use the following approximation:
 //              s = x-2.0;
 //              lgamma(x) = 0.5*s + s*P(s)/Q(s)
 //      with accuracy
 //              |P/Q - (lgamma(x)-0.5s)| < 2**-61.71
 //      Our algorithms are based on the following observation
 //
 //                             zeta(2)-1    2    zeta(3)-1    3
 // lgamma(2+s) = s*(1-Euler) + --------- * s  -  --------- * s  + ...
 //                                 2                 3
 //
 //      where Euler = 0.5772156649... is the Euler constant, which
 //      is very close to 0.5.
 //
 //   3. For x>=8, we have
 //      lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
 //      (better formula:
 //         lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
 //      Let z = 1/x, then we approximation
 //              f(z) = lgamma(x) - (x-0.5)(log(x)-1)
 //      by
 //                                  3       5             11
 //              w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
 //      where
 //              |w - f(z)| < 2**-58.74
 //
 //   4. For negative x, since (G is gamma function)
 //              -x*G(-x)*G(x) = pi/sin(pi*x),
 //      we have
 //              G(x) = pi/(sin(pi*x)*(-x)*G(-x))
 //      since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
 //      Hence, for x<0, signgam = sign(sin(pi*x)) and
 //              lgamma(x) = log(|Gamma(x)|)
 //                        = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
 //      Note: one should avoid computing pi*(-x) directly in the
 //            computation of sin(pi*(-x)).
 //
 //   5. Special Cases
 //              lgamma(2+s) ~ s*(1-Euler) for tiny s
 //              lgamma(1)=lgamma(2)=0
 //              lgamma(x) ~ -log(x) for tiny x
 //              lgamma(0) = lgamma(inf) = inf
 //              lgamma(-integer) = +-inf
 //
 //

 // Lgamma returns the natural logarithm and sign (-1 or +1) of Gamma(x).
 //
 // Special cases are:
 //	Lgamma(+Inf) = +Inf
 //	Lgamma(0) = +Inf
 //	Lgamma(-integer) = +Inf
 //	Lgamma(-Inf) = -Inf
 //	Lgamma(NaN) = NaN
 func Lgamma(x float64) (lgamma float64, sign int) {
 	const (
 		Ymin  = 1.461632144968362245
 		Two52 = 1 << 52                     // 0x4330000000000000 ~4.5036e+15
 		Two53 = 1 << 53                     // 0x4340000000000000 ~9.0072e+15
 		Two58 = 1 << 58                     // 0x4390000000000000 ~2.8823e+17
 		Tiny  = 1.0 / (1 << 70)             // 0x3b90000000000000 ~8.47033e-22
 		A0    = 7.72156649015328655494e-02  // 0x3FB3C467E37DB0C8
 		A1    = 3.22467033424113591611e-01  // 0x3FD4A34CC4A60FAD
 		A2    = 6.73523010531292681824e-02  // 0x3FB13E001A5562A7
 		A3    = 2.05808084325167332806e-02  // 0x3F951322AC92547B
 		A4    = 7.38555086081402883957e-03  // 0x3F7E404FB68FEFE8
 		A5    = 2.89051383673415629091e-03  // 0x3F67ADD8CCB7926B
 		A6    = 1.19270763183362067845e-03  // 0x3F538A94116F3F5D
 		A7    = 5.10069792153511336608e-04  // 0x3F40B6C689B99C00
 		A8    = 2.20862790713908385557e-04  // 0x3F2CF2ECED10E54D
 		A9    = 1.08011567247583939954e-04  // 0x3F1C5088987DFB07
 		A10   = 2.52144565451257326939e-05  // 0x3EFA7074428CFA52
 		A11   = 4.48640949618915160150e-05  // 0x3F07858E90A45837
 		Tc    = 1.46163214496836224576e+00  // 0x3FF762D86356BE3F
 		Tf    = -1.21486290535849611461e-01 // 0xBFBF19B9BCC38A42
 		// Tt = -(tail of Tf)
 		Tt  = -3.63867699703950536541e-18 // 0xBC50C7CAA48A971F
 		T0  = 4.83836122723810047042e-01  // 0x3FDEF72BC8EE38A2
 		T1  = -1.47587722994593911752e-01 // 0xBFC2E4278DC6C509
 		T2  = 6.46249402391333854778e-02  // 0x3FB08B4294D5419B
 		T3  = -3.27885410759859649565e-02 // 0xBFA0C9A8DF35B713
 		T4  = 1.79706750811820387126e-02  // 0x3F9266E7970AF9EC
 		T5  = -1.03142241298341437450e-02 // 0xBF851F9FBA91EC6A
 		T6  = 6.10053870246291332635e-03  // 0x3F78FCE0E370E344
 		T7  = -3.68452016781138256760e-03 // 0xBF6E2EFFB3E914D7
 		T8  = 2.25964780900612472250e-03  // 0x3F6282D32E15C915
 		T9  = -1.40346469989232843813e-03 // 0xBF56FE8EBF2D1AF1
 		T10 = 8.81081882437654011382e-04  // 0x3F4CDF0CEF61A8E9
 		T11 = -5.38595305356740546715e-04 // 0xBF41A6109C73E0EC
 		T12 = 3.15632070903625950361e-04  // 0x3F34AF6D6C0EBBF7
 		T13 = -3.12754168375120860518e-04 // 0xBF347F24ECC38C38
 		T14 = 3.35529192635519073543e-04  // 0x3F35FD3EE8C2D3F4
 		U0  = -7.72156649015328655494e-02 // 0xBFB3C467E37DB0C8
 		U1  = 6.32827064025093366517e-01  // 0x3FE4401E8B005DFF
 		U2  = 1.45492250137234768737e+00  // 0x3FF7475CD119BD6F
 		U3  = 9.77717527963372745603e-01  // 0x3FEF497644EA8450
 		U4  = 2.28963728064692451092e-01  // 0x3FCD4EAEF6010924
 		U5  = 1.33810918536787660377e-02  // 0x3F8B678BBF2BAB09
 		V1  = 2.45597793713041134822e+00  // 0x4003A5D7C2BD619C
 		V2  = 2.12848976379893395361e+00  // 0x40010725A42B18F5
 		V3  = 7.69285150456672783825e-01  // 0x3FE89DFBE45050AF
 		V4  = 1.04222645593369134254e-01  // 0x3FBAAE55D6537C88
 		V5  = 3.21709242282423911810e-03  // 0x3F6A5ABB57D0CF61
 		S0  = -7.72156649015328655494e-02 // 0xBFB3C467E37DB0C8
 		S1  = 2.14982415960608852501e-01  // 0x3FCB848B36E20878
 		S2  = 3.25778796408930981787e-01  // 0x3FD4D98F4F139F59
 		S3  = 1.46350472652464452805e-01  // 0x3FC2BB9CBEE5F2F7
 		S4  = 2.66422703033638609560e-02  // 0x3F9B481C7E939961
 		S5  = 1.84028451407337715652e-03  // 0x3F5E26B67368F239
 		S6  = 3.19475326584100867617e-05  // 0x3F00BFECDD17E945
 		R1  = 1.39200533467621045958e+00  // 0x3FF645A762C4AB74
 		R2  = 7.21935547567138069525e-01  // 0x3FE71A1893D3DCDC
 		R3  = 1.71933865632803078993e-01  // 0x3FC601EDCCFBDF27
 		R4  = 1.86459191715652901344e-02  // 0x3F9317EA742ED475
 		R5  = 7.77942496381893596434e-04  // 0x3F497DDACA41A95B
 		R6  = 7.32668430744625636189e-06  // 0x3EDEBAF7A5B38140
 		W0  = 4.18938533204672725052e-01  // 0x3FDACFE390C97D69
 		W1  = 8.33333333333329678849e-02  // 0x3FB555555555553B
 		W2  = -2.77777777728775536470e-03 // 0xBF66C16C16B02E5C
 		W3  = 7.93650558643019558500e-04  // 0x3F4A019F98CF38B6
 		W4  = -5.95187557450339963135e-04 // 0xBF4380CB8C0FE741
 		W5  = 8.36339918996282139126e-04  // 0x3F4B67BA4CDAD5D1
 		W6  = -1.63092934096575273989e-03 // 0xBF5AB89D0B9E43E4
 	)
 	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
 	// when compiler does it for us
 	// special cases
 	sign = 1
 	switch {
 	case x != x: // IsNaN(x):
 		lgamma = x
 		return
 	case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0):
 		lgamma = x
 		return
 	case x == 0:
 		lgamma = Inf(1)
 		return
 	}

 	neg := false
 	if x < 0 {
 		x = -x
 		neg = true
 	}

 	if x < Tiny { // if |x| < 2**-70, return -log(|x|)
 		if neg {
 			sign = -1
 		}
 		lgamma = -Log(x)
 		return
 	}
 	var nadj float64
 	if neg {
 		if x >= Two52 { // |x| >= 2**52, must be -integer
 			lgamma = Inf(1)
 			return
 		}
 		t := sinPi(x)
 		if t == 0 {
 			lgamma = Inf(1) // -integer
 			return
 		}
 		nadj = Log(Pi / Fabs(t*x))
 		if t < 0 {
 			sign = -1
 		}
 	}

 	switch {
 	case x == 1 || x == 2: // purge off 1 and 2
 		lgamma = 0
 		return
 	case x < 2: // use lgamma(x) = lgamma(x+1) - log(x)
 		var y float64
 		var i int
 		if x <= 0.9 {
 			lgamma = -Log(x)
 			switch {
 			case x >= (Ymin - 1 + 0.27): // 0.7316 <= x <=  0.9
 				y = 1 - x
 				i = 0
 			case x >= (Ymin - 1 - 0.27): // 0.2316 <= x < 0.7316
 				y = x - (Tc - 1)
 				i = 1
 			default: // 0 < x < 0.2316
 				y = x
 				i = 2
 			}
 		} else {
 			lgamma = 0
 			switch {
 			case x >= (Ymin + 0.27): // 1.7316 <= x < 2
 				y = 2 - x
 				i = 0
 			case x >= (Ymin - 0.27): // 1.2316 <= x < 1.7316
 				y = x - Tc
 				i = 1
 			default: // 0.9 < x < 1.2316
 				y = x - 1
 				i = 2
 			}
 		}
 		switch i {
 		case 0:
 			z := y * y
 			p1 := A0 + z*(A2+z*(A4+z*(A6+z*(A8+z*A10))))
 			p2 := z * (A1 + z*(A3+z*(A5+z*(A7+z*(A9+z*A11)))))
 			p := y*p1 + p2
 			lgamma += (p - 0.5*y)
 		case 1:
 			z := y * y
 			w := z * y
 			p1 := T0 + w*(T3+w*(T6+w*(T9+w*T12))) // parallel comp
 			p2 := T1 + w*(T4+w*(T7+w*(T10+w*T13)))
 			p3 := T2 + w*(T5+w*(T8+w*(T11+w*T14)))
 			p := z*p1 - (Tt - w*(p2+y*p3))
 			lgamma += (Tf + p)
 		case 2:
 			p1 := y * (U0 + y*(U1+y*(U2+y*(U3+y*(U4+y*U5)))))
 			p2 := 1 + y*(V1+y*(V2+y*(V3+y*(V4+y*V5))))
 			lgamma += (-0.5*y + p1/p2)
 		}
 	case x < 8: // 2 <= x < 8
 		i := int(x)
 		y := x - float64(i)
 		p := y * (S0 + y*(S1+y*(S2+y*(S3+y*(S4+y*(S5+y*S6))))))
 		q := 1 + y*(R1+y*(R2+y*(R3+y*(R4+y*(R5+y*R6)))))
 		lgamma = 0.5*y + p/q
 		z := 1.0 // Lgamma(1+s) = Log(s) + Lgamma(s)
 		switch i {
 		case 7:
 			z *= (y + 6)
 			fallthrough
 		case 6:
 			z *= (y + 5)
 			fallthrough
 		case 5:
 			z *= (y + 4)
 			fallthrough
 		case 4:
 			z *= (y + 3)
 			fallthrough
 		case 3:
 			z *= (y + 2)
 			lgamma += Log(z)
 		}
 	case x < Two58: // 8 <= x < 2**58
 		t := Log(x)
 		z := 1 / x
 		y := z * z
 		w := W0 + z*(W1+y*(W2+y*(W3+y*(W4+y*(W5+y*W6)))))
 		lgamma = (x-0.5)*(t-1) + w
 	default: // 2**58 <= x <= Inf
 		lgamma = x * (Log(x) - 1)
 	}
 	if neg {
 		lgamma = nadj - lgamma
 	}
 	return
 }

 // sinPi(x) is a helper function for negative x
 func sinPi(x float64) float64 {
 	const (
 		Two52 = 1 << 52 // 0x4330000000000000 ~4.5036e+15
 		Two53 = 1 << 53 // 0x4340000000000000 ~9.0072e+15
 	)
 	if x < 0.25 {
 		return -Sin(Pi * x)
 	}

 	// argument reduction
 	z := Floor(x)
 	var n int
 	if z != x { // inexact
 		x = Fmod(x, 2)
 		n = int(x * 4)
 	} else {
 		if x >= Two53 { // x must be even
 			x = 0
 			n = 0
 		} else {
 			if x < Two52 {
 				z = x + Two52 // exact
 			}
 			n = int(1 & Float64bits(z))
 			x = float64(n)
 			n <<= 2
 		}
 	}
 	switch n {
 	case 0:
 		x = Sin(Pi * x)
 	case 1, 2:
 		x = Cos(Pi * (0.5 - x))
 	case 3, 4:
 		x = Sin(Pi * (1 - x))
 	case 5, 6:
 		x = -Cos(Pi * (x - 1.5))
 	default:
 		x = Sin(Pi * (x - 2))
 	}
 	return -x
 }
	// 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

	/*
	Floating-point logarithm of the Gamma function.
	*/

	// The original C code and the long comment below are
	// from FreeBSD's /usr/src/lib/msun/src/e_lgamma_r.c and
	// came with this notice. The go code is a simplified
	// version of the original C.
	//
	// ====================================================
	// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
	//
	// Developed at SunPro, a Sun Microsystems, Inc. business.
	// Permission to use, copy, modify, and distribute this
	// software is freely granted, provided that this notice
	// is preserved.
	// ====================================================
	//
	// __ieee754_lgamma_r(x, signgamp)
	// Reentrant version of the logarithm of the Gamma function
	// with user provided pointer for the sign of Gamma(x).
	//
	// Method:
	// 1. Argument Reduction for 0 < x <= 8
	// Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
	// reduce x to a number in [1.5,2.5] by
	// lgamma(1+s) = log(s) + lgamma(s)
	// for example,
	// lgamma(7.3) = log(6.3) + lgamma(6.3)
	// = log(6.3*5.3) + lgamma(5.3)
	// = log(6.35.34.33.32.3) + lgamma(2.3)
	// 2. Polynomial approximation of lgamma around its
	// minimum (ymin=1.461632144968362245) to maintain monotonicity.
	// On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
	// Let z = x-ymin;
	// lgamma(x) = -1.214862905358496078218 + z*2poly(z)
	// poly(z) is a 14 degree polynomial.
	// 2. Rational approximation in the primary interval [2,3]
	// We use the following approximation:
	// s = x-2.0;
	// lgamma(x) = 0.5s + sP(s)/Q(s)
	// with accuracy
	// \|P/Q - (lgamma(x)-0.5s)\| < 2**-61.71
	// Our algorithms are based on the following observation
	//
	// zeta(2)-1 2 zeta(3)-1 3
	// lgamma(2+s) = s(1-Euler) + --------- s - --------- * s + ...
	// 2 3
	//
	// where Euler = 0.5772156649... is the Euler constant, which
	// is very close to 0.5.
	//
	// 3. For x>=8, we have
	// lgamma(x)~(x-0.5)log(x)-x+0.5log(2pi)+1/(12x)-1/(360x*3)+....
	// (better formula:
	// lgamma(x)~(x-0.5)(log(x)-1)-.5(log(2pi)-1) + ...)
	// Let z = 1/x, then we approximation
	// f(z) = lgamma(x) - (x-0.5)(log(x)-1)
	// by
	// 3 5 11
	// w = w0 + w1z + w2z + w3z + ... + w6z
	// where
	// \|w - f(z)\| < 2**-58.74
	//
	// 4. For negative x, since (G is gamma function)
	// -xG(-x)G(x) = pi/sin(pi*x),
	// we have
	// G(x) = pi/(sin(pix)(-x)*G(-x))
	// since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
	// Hence, for x<0, signgam = sign(sin(pi*x)) and
	// lgamma(x) = log(\|Gamma(x)\|)
	// = log(pi/(\|xsin(pix)\|)) - lgamma(-x);
	// Note: one should avoid computing pi*(-x) directly in the
	// computation of sin(pi*(-x)).
	//
	// 5. Special Cases
	// lgamma(2+s) ~ s*(1-Euler) for tiny s
	// lgamma(1)=lgamma(2)=0
	// lgamma(x) ~ -log(x) for tiny x
	// lgamma(0) = lgamma(inf) = inf
	// lgamma(-integer) = +-inf
	//
	//

	// Lgamma returns the natural logarithm and sign (-1 or +1) of Gamma(x).
	//
	// Special cases are:
	// Lgamma(+Inf) = +Inf
	// Lgamma(0) = +Inf
	// Lgamma(-integer) = +Inf
	// Lgamma(-Inf) = -Inf
	// Lgamma(NaN) = NaN
	func Lgamma(x float64) (lgamma float64, sign int) {
	const (
	Ymin = 1.461632144968362245
	Two52 = 1 << 52 // 0x4330000000000000 ~4.5036e+15
	Two53 = 1 << 53 // 0x4340000000000000 ~9.0072e+15
	Two58 = 1 << 58 // 0x4390000000000000 ~2.8823e+17
	Tiny = 1.0 / (1 << 70) // 0x3b90000000000000 ~8.47033e-22
	A0 = 7.72156649015328655494e-02 // 0x3FB3C467E37DB0C8
	A1 = 3.22467033424113591611e-01 // 0x3FD4A34CC4A60FAD
	A2 = 6.73523010531292681824e-02 // 0x3FB13E001A5562A7
	A3 = 2.05808084325167332806e-02 // 0x3F951322AC92547B
	A4 = 7.38555086081402883957e-03 // 0x3F7E404FB68FEFE8
	A5 = 2.89051383673415629091e-03 // 0x3F67ADD8CCB7926B
	A6 = 1.19270763183362067845e-03 // 0x3F538A94116F3F5D
	A7 = 5.10069792153511336608e-04 // 0x3F40B6C689B99C00
	A8 = 2.20862790713908385557e-04 // 0x3F2CF2ECED10E54D
	A9 = 1.08011567247583939954e-04 // 0x3F1C5088987DFB07
	A10 = 2.52144565451257326939e-05 // 0x3EFA7074428CFA52
	A11 = 4.48640949618915160150e-05 // 0x3F07858E90A45837
	Tc = 1.46163214496836224576e+00 // 0x3FF762D86356BE3F
	Tf = -1.21486290535849611461e-01 // 0xBFBF19B9BCC38A42
	// Tt = -(tail of Tf)
	Tt = -3.63867699703950536541e-18 // 0xBC50C7CAA48A971F
	T0 = 4.83836122723810047042e-01 // 0x3FDEF72BC8EE38A2
	T1 = -1.47587722994593911752e-01 // 0xBFC2E4278DC6C509
	T2 = 6.46249402391333854778e-02 // 0x3FB08B4294D5419B
	T3 = -3.27885410759859649565e-02 // 0xBFA0C9A8DF35B713
	T4 = 1.79706750811820387126e-02 // 0x3F9266E7970AF9EC
	T5 = -1.03142241298341437450e-02 // 0xBF851F9FBA91EC6A
	T6 = 6.10053870246291332635e-03 // 0x3F78FCE0E370E344
	T7 = -3.68452016781138256760e-03 // 0xBF6E2EFFB3E914D7
	T8 = 2.25964780900612472250e-03 // 0x3F6282D32E15C915
	T9 = -1.40346469989232843813e-03 // 0xBF56FE8EBF2D1AF1
	T10 = 8.81081882437654011382e-04 // 0x3F4CDF0CEF61A8E9
	T11 = -5.38595305356740546715e-04 // 0xBF41A6109C73E0EC
	T12 = 3.15632070903625950361e-04 // 0x3F34AF6D6C0EBBF7
	T13 = -3.12754168375120860518e-04 // 0xBF347F24ECC38C38
	T14 = 3.35529192635519073543e-04 // 0x3F35FD3EE8C2D3F4
	U0 = -7.72156649015328655494e-02 // 0xBFB3C467E37DB0C8
	U1 = 6.32827064025093366517e-01 // 0x3FE4401E8B005DFF
	U2 = 1.45492250137234768737e+00 // 0x3FF7475CD119BD6F
	U3 = 9.77717527963372745603e-01 // 0x3FEF497644EA8450
	U4 = 2.28963728064692451092e-01 // 0x3FCD4EAEF6010924
	U5 = 1.33810918536787660377e-02 // 0x3F8B678BBF2BAB09
	V1 = 2.45597793713041134822e+00 // 0x4003A5D7C2BD619C
	V2 = 2.12848976379893395361e+00 // 0x40010725A42B18F5
	V3 = 7.69285150456672783825e-01 // 0x3FE89DFBE45050AF
	V4 = 1.04222645593369134254e-01 // 0x3FBAAE55D6537C88
	V5 = 3.21709242282423911810e-03 // 0x3F6A5ABB57D0CF61
	S0 = -7.72156649015328655494e-02 // 0xBFB3C467E37DB0C8
	S1 = 2.14982415960608852501e-01 // 0x3FCB848B36E20878
	S2 = 3.25778796408930981787e-01 // 0x3FD4D98F4F139F59
	S3 = 1.46350472652464452805e-01 // 0x3FC2BB9CBEE5F2F7
	S4 = 2.66422703033638609560e-02 // 0x3F9B481C7E939961
	S5 = 1.84028451407337715652e-03 // 0x3F5E26B67368F239
	S6 = 3.19475326584100867617e-05 // 0x3F00BFECDD17E945
	R1 = 1.39200533467621045958e+00 // 0x3FF645A762C4AB74
	R2 = 7.21935547567138069525e-01 // 0x3FE71A1893D3DCDC
	R3 = 1.71933865632803078993e-01 // 0x3FC601EDCCFBDF27
	R4 = 1.86459191715652901344e-02 // 0x3F9317EA742ED475
	R5 = 7.77942496381893596434e-04 // 0x3F497DDACA41A95B
	R6 = 7.32668430744625636189e-06 // 0x3EDEBAF7A5B38140
	W0 = 4.18938533204672725052e-01 // 0x3FDACFE390C97D69
	W1 = 8.33333333333329678849e-02 // 0x3FB555555555553B
	W2 = -2.77777777728775536470e-03 // 0xBF66C16C16B02E5C
	W3 = 7.93650558643019558500e-04 // 0x3F4A019F98CF38B6
	W4 = -5.95187557450339963135e-04 // 0xBF4380CB8C0FE741
	W5 = 8.36339918996282139126e-04 // 0x3F4B67BA4CDAD5D1
	W6 = -1.63092934096575273989e-03 // 0xBF5AB89D0B9E43E4
	)
	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
	// when compiler does it for us
	// special cases
	sign = 1
	switch {
	case x != x: // IsNaN(x):
	lgamma = x
	return
	case x < -MaxFloat64 \|\| x > MaxFloat64: // IsInf(x, 0):
	lgamma = x
	return
	case x == 0:
	lgamma = Inf(1)
	return
	}

	neg := false
	if x < 0 {
	x = -x
	neg = true
	}

	if x < Tiny { // if \|x\| < 2**-70, return -log(\|x\|)
	if neg {
	sign = -1
	}
	lgamma = -Log(x)
	return
	}
	var nadj float64
	if neg {
	if x >= Two52 { // \|x\| >= 2**52, must be -integer
	lgamma = Inf(1)
	return
	}
	t := sinPi(x)
	if t == 0 {
	lgamma = Inf(1) // -integer
	return
	}
	nadj = Log(Pi / Fabs(t*x))
	if t < 0 {
	sign = -1
	}
	}

	switch {
	case x == 1 \|\| x == 2: // purge off 1 and 2
	lgamma = 0
	return
	case x < 2: // use lgamma(x) = lgamma(x+1) - log(x)
	var y float64
	var i int
	if x <= 0.9 {
	lgamma = -Log(x)
	switch {
	case x >= (Ymin - 1 + 0.27): // 0.7316 <= x <= 0.9
	y = 1 - x
	i = 0
	case x >= (Ymin - 1 - 0.27): // 0.2316 <= x < 0.7316
	y = x - (Tc - 1)
	i = 1
	default: // 0 < x < 0.2316
	y = x
	i = 2
	}
	} else {
	lgamma = 0
	switch {
	case x >= (Ymin + 0.27): // 1.7316 <= x < 2
	y = 2 - x
	i = 0
	case x >= (Ymin - 0.27): // 1.2316 <= x < 1.7316
	y = x - Tc
	i = 1
	default: // 0.9 < x < 1.2316
	y = x - 1
	i = 2
	}
	}
	switch i {
	case 0:
	z := y * y
	p1 := A0 + z(A2+z(A4+z(A6+z(A8+z*A10))))
	p2 := z * (A1 + z(A3+z(A5+z(A7+z(A9+z*A11)))))
	p := y*p1 + p2
	lgamma += (p - 0.5*y)
	case 1:
	z := y * y
	w := z * y
	p1 := T0 + w(T3+w(T6+w(T9+wT12))) // parallel comp
	p2 := T1 + w(T4+w(T7+w(T10+wT13)))
	p3 := T2 + w(T5+w(T8+w(T11+wT14)))
	p := zp1 - (Tt - w(p2+y*p3))
	lgamma += (Tf + p)
	case 2:
	p1 := y * (U0 + y(U1+y(U2+y(U3+y(U4+y*U5)))))
	p2 := 1 + y(V1+y(V2+y(V3+y(V4+y*V5))))
	lgamma += (-0.5*y + p1/p2)
	}
	case x < 8: // 2 <= x < 8
	i := int(x)
	y := x - float64(i)
	p := y * (S0 + y(S1+y(S2+y(S3+y(S4+y(S5+yS6))))))
	q := 1 + y(R1+y(R2+y(R3+y(R4+y(R5+yR6)))))
	lgamma = 0.5*y + p/q
	z := 1.0 // Lgamma(1+s) = Log(s) + Lgamma(s)
	switch i {
	case 7:
	z *= (y + 6)
	fallthrough
	case 6:
	z *= (y + 5)
	fallthrough
	case 5:
	z *= (y + 4)
	fallthrough
	case 4:
	z *= (y + 3)
	fallthrough
	case 3:
	z *= (y + 2)
	lgamma += Log(z)
	}
	case x < Two58: // 8 <= x < 2**58
	t := Log(x)
	z := 1 / x
	y := z * z
	w := W0 + z(W1+y(W2+y(W3+y(W4+y(W5+yW6)))))
	lgamma = (x-0.5)*(t-1) + w
	default: // 2**58 <= x <= Inf
	lgamma = x * (Log(x) - 1)
	}
	if neg {
	lgamma = nadj - lgamma
	}
	return
	}

	// sinPi(x) is a helper function for negative x
	func sinPi(x float64) float64 {
	const (
	Two52 = 1 << 52 // 0x4330000000000000 ~4.5036e+15
	Two53 = 1 << 53 // 0x4340000000000000 ~9.0072e+15
	)
	if x < 0.25 {
	return -Sin(Pi * x)
	}

	// argument reduction
	z := Floor(x)
	var n int
	if z != x { // inexact
	x = Fmod(x, 2)
	n = int(x * 4)
	} else {
	if x >= Two53 { // x must be even
	x = 0
	n = 0
	} else {
	if x < Two52 {
	z = x + Two52 // exact
	}
	n = int(1 & Float64bits(z))
	x = float64(n)
	n <<= 2
	}
	}
	switch n {
	case 0:
	x = Sin(Pi * x)
	case 1, 2:
	x = Cos(Pi * (0.5 - x))
	case 3, 4:
	x = Sin(Pi * (1 - x))
	case 5, 6:
	x = -Cos(Pi * (x - 1.5))
	default:
	x = Sin(Pi * (x - 2))
	}
	return -x
	}