Source file src/math/lgamma.go

Documentation: math

		 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  /*
		 8  	Floating-point logarithm of the Gamma function.
		 9  */
		10  
		11  // The original C code and the long comment below are
		12  // from FreeBSD's /usr/src/lib/msun/src/e_lgamma_r.c and
		13  // came with this notice. The go code is a simplified
		14  // version of the original C.
		15  //
		16  // ====================================================
		17  // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
		18  //
		19  // Developed at SunPro, a Sun Microsystems, Inc. business.
		20  // Permission to use, copy, modify, and distribute this
		21  // software is freely granted, provided that this notice
		22  // is preserved.
		23  // ====================================================
		24  //
		25  // __ieee754_lgamma_r(x, signgamp)
		26  // Reentrant version of the logarithm of the Gamma function
		27  // with user provided pointer for the sign of Gamma(x).
		28  //
		29  // Method:
		30  //	 1. Argument Reduction for 0 < x <= 8
		31  //			Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
		32  //			reduce x to a number in [1.5,2.5] by
		33  //							lgamma(1+s) = log(s) + lgamma(s)
		34  //			for example,
		35  //							lgamma(7.3) = log(6.3) + lgamma(6.3)
		36  //													= log(6.3*5.3) + lgamma(5.3)
		37  //													= log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
		38  //	 2. Polynomial approximation of lgamma around its
		39  //			minimum (ymin=1.461632144968362245) to maintain monotonicity.
		40  //			On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
		41  //							Let z = x-ymin;
		42  //							lgamma(x) = -1.214862905358496078218 + z**2*poly(z)
		43  //							poly(z) is a 14 degree polynomial.
		44  //	 2. Rational approximation in the primary interval [2,3]
		45  //			We use the following approximation:
		46  //							s = x-2.0;
		47  //							lgamma(x) = 0.5*s + s*P(s)/Q(s)
		48  //			with accuracy
		49  //							|P/Q - (lgamma(x)-0.5s)| < 2**-61.71
		50  //			Our algorithms are based on the following observation
		51  //
		52  //														 zeta(2)-1		2		zeta(3)-1		3
		53  // lgamma(2+s) = s*(1-Euler) + --------- * s	-	--------- * s	+ ...
		54  //																 2								 3
		55  //
		56  //			where Euler = 0.5772156649... is the Euler constant, which
		57  //			is very close to 0.5.
		58  //
		59  //	 3. For x>=8, we have
		60  //			lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
		61  //			(better formula:
		62  //				 lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
		63  //			Let z = 1/x, then we approximation
		64  //							f(z) = lgamma(x) - (x-0.5)(log(x)-1)
		65  //			by
		66  //																	3			 5						 11
		67  //							w = w0 + w1*z + w2*z	+ w3*z	+ ... + w6*z
		68  //			where
		69  //							|w - f(z)| < 2**-58.74
		70  //
		71  //	 4. For negative x, since (G is gamma function)
		72  //							-x*G(-x)*G(x) = pi/sin(pi*x),
		73  //			we have
		74  //							G(x) = pi/(sin(pi*x)*(-x)*G(-x))
		75  //			since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
		76  //			Hence, for x<0, signgam = sign(sin(pi*x)) and
		77  //							lgamma(x) = log(|Gamma(x)|)
		78  //												= log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
		79  //			Note: one should avoid computing pi*(-x) directly in the
		80  //						computation of sin(pi*(-x)).
		81  //
		82  //	 5. Special Cases
		83  //							lgamma(2+s) ~ s*(1-Euler) for tiny s
		84  //							lgamma(1)=lgamma(2)=0
		85  //							lgamma(x) ~ -log(x) for tiny x
		86  //							lgamma(0) = lgamma(inf) = inf
		87  //							lgamma(-integer) = +-inf
		88  //
		89  //
		90  
		91  var _lgamA = [...]float64{
		92  	7.72156649015328655494e-02, // 0x3FB3C467E37DB0C8
		93  	3.22467033424113591611e-01, // 0x3FD4A34CC4A60FAD
		94  	6.73523010531292681824e-02, // 0x3FB13E001A5562A7
		95  	2.05808084325167332806e-02, // 0x3F951322AC92547B
		96  	7.38555086081402883957e-03, // 0x3F7E404FB68FEFE8
		97  	2.89051383673415629091e-03, // 0x3F67ADD8CCB7926B
		98  	1.19270763183362067845e-03, // 0x3F538A94116F3F5D
		99  	5.10069792153511336608e-04, // 0x3F40B6C689B99C00
	 100  	2.20862790713908385557e-04, // 0x3F2CF2ECED10E54D
	 101  	1.08011567247583939954e-04, // 0x3F1C5088987DFB07
	 102  	2.52144565451257326939e-05, // 0x3EFA7074428CFA52
	 103  	4.48640949618915160150e-05, // 0x3F07858E90A45837
	 104  }
	 105  var _lgamR = [...]float64{
	 106  	1.0,												// placeholder
	 107  	1.39200533467621045958e+00, // 0x3FF645A762C4AB74
	 108  	7.21935547567138069525e-01, // 0x3FE71A1893D3DCDC
	 109  	1.71933865632803078993e-01, // 0x3FC601EDCCFBDF27
	 110  	1.86459191715652901344e-02, // 0x3F9317EA742ED475
	 111  	7.77942496381893596434e-04, // 0x3F497DDACA41A95B
	 112  	7.32668430744625636189e-06, // 0x3EDEBAF7A5B38140
	 113  }
	 114  var _lgamS = [...]float64{
	 115  	-7.72156649015328655494e-02, // 0xBFB3C467E37DB0C8
	 116  	2.14982415960608852501e-01,	// 0x3FCB848B36E20878
	 117  	3.25778796408930981787e-01,	// 0x3FD4D98F4F139F59
	 118  	1.46350472652464452805e-01,	// 0x3FC2BB9CBEE5F2F7
	 119  	2.66422703033638609560e-02,	// 0x3F9B481C7E939961
	 120  	1.84028451407337715652e-03,	// 0x3F5E26B67368F239
	 121  	3.19475326584100867617e-05,	// 0x3F00BFECDD17E945
	 122  }
	 123  var _lgamT = [...]float64{
	 124  	4.83836122723810047042e-01,	// 0x3FDEF72BC8EE38A2
	 125  	-1.47587722994593911752e-01, // 0xBFC2E4278DC6C509
	 126  	6.46249402391333854778e-02,	// 0x3FB08B4294D5419B
	 127  	-3.27885410759859649565e-02, // 0xBFA0C9A8DF35B713
	 128  	1.79706750811820387126e-02,	// 0x3F9266E7970AF9EC
	 129  	-1.03142241298341437450e-02, // 0xBF851F9FBA91EC6A
	 130  	6.10053870246291332635e-03,	// 0x3F78FCE0E370E344
	 131  	-3.68452016781138256760e-03, // 0xBF6E2EFFB3E914D7
	 132  	2.25964780900612472250e-03,	// 0x3F6282D32E15C915
	 133  	-1.40346469989232843813e-03, // 0xBF56FE8EBF2D1AF1
	 134  	8.81081882437654011382e-04,	// 0x3F4CDF0CEF61A8E9
	 135  	-5.38595305356740546715e-04, // 0xBF41A6109C73E0EC
	 136  	3.15632070903625950361e-04,	// 0x3F34AF6D6C0EBBF7
	 137  	-3.12754168375120860518e-04, // 0xBF347F24ECC38C38
	 138  	3.35529192635519073543e-04,	// 0x3F35FD3EE8C2D3F4
	 139  }
	 140  var _lgamU = [...]float64{
	 141  	-7.72156649015328655494e-02, // 0xBFB3C467E37DB0C8
	 142  	6.32827064025093366517e-01,	// 0x3FE4401E8B005DFF
	 143  	1.45492250137234768737e+00,	// 0x3FF7475CD119BD6F
	 144  	9.77717527963372745603e-01,	// 0x3FEF497644EA8450
	 145  	2.28963728064692451092e-01,	// 0x3FCD4EAEF6010924
	 146  	1.33810918536787660377e-02,	// 0x3F8B678BBF2BAB09
	 147  }
	 148  var _lgamV = [...]float64{
	 149  	1.0,
	 150  	2.45597793713041134822e+00, // 0x4003A5D7C2BD619C
	 151  	2.12848976379893395361e+00, // 0x40010725A42B18F5
	 152  	7.69285150456672783825e-01, // 0x3FE89DFBE45050AF
	 153  	1.04222645593369134254e-01, // 0x3FBAAE55D6537C88
	 154  	3.21709242282423911810e-03, // 0x3F6A5ABB57D0CF61
	 155  }
	 156  var _lgamW = [...]float64{
	 157  	4.18938533204672725052e-01,	// 0x3FDACFE390C97D69
	 158  	8.33333333333329678849e-02,	// 0x3FB555555555553B
	 159  	-2.77777777728775536470e-03, // 0xBF66C16C16B02E5C
	 160  	7.93650558643019558500e-04,	// 0x3F4A019F98CF38B6
	 161  	-5.95187557450339963135e-04, // 0xBF4380CB8C0FE741
	 162  	8.36339918996282139126e-04,	// 0x3F4B67BA4CDAD5D1
	 163  	-1.63092934096575273989e-03, // 0xBF5AB89D0B9E43E4
	 164  }
	 165  
	 166  // Lgamma returns the natural logarithm and sign (-1 or +1) of Gamma(x).
	 167  //
	 168  // Special cases are:
	 169  //	Lgamma(+Inf) = +Inf
	 170  //	Lgamma(0) = +Inf
	 171  //	Lgamma(-integer) = +Inf
	 172  //	Lgamma(-Inf) = -Inf
	 173  //	Lgamma(NaN) = NaN
	 174  func Lgamma(x float64) (lgamma float64, sign int) {
	 175  	const (
	 176  		Ymin	= 1.461632144968362245
	 177  		Two52 = 1 << 52										 // 0x4330000000000000 ~4.5036e+15
	 178  		Two53 = 1 << 53										 // 0x4340000000000000 ~9.0072e+15
	 179  		Two58 = 1 << 58										 // 0x4390000000000000 ~2.8823e+17
	 180  		Tiny	= 1.0 / (1 << 70)						 // 0x3b90000000000000 ~8.47033e-22
	 181  		Tc		= 1.46163214496836224576e+00	// 0x3FF762D86356BE3F
	 182  		Tf		= -1.21486290535849611461e-01 // 0xBFBF19B9BCC38A42
	 183  		// Tt = -(tail of Tf)
	 184  		Tt = -3.63867699703950536541e-18 // 0xBC50C7CAA48A971F
	 185  	)
	 186  	// special cases
	 187  	sign = 1
	 188  	switch {
	 189  	case IsNaN(x):
	 190  		lgamma = x
	 191  		return
	 192  	case IsInf(x, 0):
	 193  		lgamma = x
	 194  		return
	 195  	case x == 0:
	 196  		lgamma = Inf(1)
	 197  		return
	 198  	}
	 199  
	 200  	neg := false
	 201  	if x < 0 {
	 202  		x = -x
	 203  		neg = true
	 204  	}
	 205  
	 206  	if x < Tiny { // if |x| < 2**-70, return -log(|x|)
	 207  		if neg {
	 208  			sign = -1
	 209  		}
	 210  		lgamma = -Log(x)
	 211  		return
	 212  	}
	 213  	var nadj float64
	 214  	if neg {
	 215  		if x >= Two52 { // |x| >= 2**52, must be -integer
	 216  			lgamma = Inf(1)
	 217  			return
	 218  		}
	 219  		t := sinPi(x)
	 220  		if t == 0 {
	 221  			lgamma = Inf(1) // -integer
	 222  			return
	 223  		}
	 224  		nadj = Log(Pi / Abs(t*x))
	 225  		if t < 0 {
	 226  			sign = -1
	 227  		}
	 228  	}
	 229  
	 230  	switch {
	 231  	case x == 1 || x == 2: // purge off 1 and 2
	 232  		lgamma = 0
	 233  		return
	 234  	case x < 2: // use lgamma(x) = lgamma(x+1) - log(x)
	 235  		var y float64
	 236  		var i int
	 237  		if x <= 0.9 {
	 238  			lgamma = -Log(x)
	 239  			switch {
	 240  			case x >= (Ymin - 1 + 0.27): // 0.7316 <= x <=	0.9
	 241  				y = 1 - x
	 242  				i = 0
	 243  			case x >= (Ymin - 1 - 0.27): // 0.2316 <= x < 0.7316
	 244  				y = x - (Tc - 1)
	 245  				i = 1
	 246  			default: // 0 < x < 0.2316
	 247  				y = x
	 248  				i = 2
	 249  			}
	 250  		} else {
	 251  			lgamma = 0
	 252  			switch {
	 253  			case x >= (Ymin + 0.27): // 1.7316 <= x < 2
	 254  				y = 2 - x
	 255  				i = 0
	 256  			case x >= (Ymin - 0.27): // 1.2316 <= x < 1.7316
	 257  				y = x - Tc
	 258  				i = 1
	 259  			default: // 0.9 < x < 1.2316
	 260  				y = x - 1
	 261  				i = 2
	 262  			}
	 263  		}
	 264  		switch i {
	 265  		case 0:
	 266  			z := y * y
	 267  			p1 := _lgamA[0] + z*(_lgamA[2]+z*(_lgamA[4]+z*(_lgamA[6]+z*(_lgamA[8]+z*_lgamA[10]))))
	 268  			p2 := z * (_lgamA[1] + z*(+_lgamA[3]+z*(_lgamA[5]+z*(_lgamA[7]+z*(_lgamA[9]+z*_lgamA[11])))))
	 269  			p := y*p1 + p2
	 270  			lgamma += (p - 0.5*y)
	 271  		case 1:
	 272  			z := y * y
	 273  			w := z * y
	 274  			p1 := _lgamT[0] + w*(_lgamT[3]+w*(_lgamT[6]+w*(_lgamT[9]+w*_lgamT[12]))) // parallel comp
	 275  			p2 := _lgamT[1] + w*(_lgamT[4]+w*(_lgamT[7]+w*(_lgamT[10]+w*_lgamT[13])))
	 276  			p3 := _lgamT[2] + w*(_lgamT[5]+w*(_lgamT[8]+w*(_lgamT[11]+w*_lgamT[14])))
	 277  			p := z*p1 - (Tt - w*(p2+y*p3))
	 278  			lgamma += (Tf + p)
	 279  		case 2:
	 280  			p1 := y * (_lgamU[0] + y*(_lgamU[1]+y*(_lgamU[2]+y*(_lgamU[3]+y*(_lgamU[4]+y*_lgamU[5])))))
	 281  			p2 := 1 + y*(_lgamV[1]+y*(_lgamV[2]+y*(_lgamV[3]+y*(_lgamV[4]+y*_lgamV[5]))))
	 282  			lgamma += (-0.5*y + p1/p2)
	 283  		}
	 284  	case x < 8: // 2 <= x < 8
	 285  		i := int(x)
	 286  		y := x - float64(i)
	 287  		p := y * (_lgamS[0] + y*(_lgamS[1]+y*(_lgamS[2]+y*(_lgamS[3]+y*(_lgamS[4]+y*(_lgamS[5]+y*_lgamS[6]))))))
	 288  		q := 1 + y*(_lgamR[1]+y*(_lgamR[2]+y*(_lgamR[3]+y*(_lgamR[4]+y*(_lgamR[5]+y*_lgamR[6])))))
	 289  		lgamma = 0.5*y + p/q
	 290  		z := 1.0 // Lgamma(1+s) = Log(s) + Lgamma(s)
	 291  		switch i {
	 292  		case 7:
	 293  			z *= (y + 6)
	 294  			fallthrough
	 295  		case 6:
	 296  			z *= (y + 5)
	 297  			fallthrough
	 298  		case 5:
	 299  			z *= (y + 4)
	 300  			fallthrough
	 301  		case 4:
	 302  			z *= (y + 3)
	 303  			fallthrough
	 304  		case 3:
	 305  			z *= (y + 2)
	 306  			lgamma += Log(z)
	 307  		}
	 308  	case x < Two58: // 8 <= x < 2**58
	 309  		t := Log(x)
	 310  		z := 1 / x
	 311  		y := z * z
	 312  		w := _lgamW[0] + z*(_lgamW[1]+y*(_lgamW[2]+y*(_lgamW[3]+y*(_lgamW[4]+y*(_lgamW[5]+y*_lgamW[6])))))
	 313  		lgamma = (x-0.5)*(t-1) + w
	 314  	default: // 2**58 <= x <= Inf
	 315  		lgamma = x * (Log(x) - 1)
	 316  	}
	 317  	if neg {
	 318  		lgamma = nadj - lgamma
	 319  	}
	 320  	return
	 321  }
	 322  
	 323  // sinPi(x) is a helper function for negative x
	 324  func sinPi(x float64) float64 {
	 325  	const (
	 326  		Two52 = 1 << 52 // 0x4330000000000000 ~4.5036e+15
	 327  		Two53 = 1 << 53 // 0x4340000000000000 ~9.0072e+15
	 328  	)
	 329  	if x < 0.25 {
	 330  		return -Sin(Pi * x)
	 331  	}
	 332  
	 333  	// argument reduction
	 334  	z := Floor(x)
	 335  	var n int
	 336  	if z != x { // inexact
	 337  		x = Mod(x, 2)
	 338  		n = int(x * 4)
	 339  	} else {
	 340  		if x >= Two53 { // x must be even
	 341  			x = 0
	 342  			n = 0
	 343  		} else {
	 344  			if x < Two52 {
	 345  				z = x + Two52 // exact
	 346  			}
	 347  			n = int(1 & Float64bits(z))
	 348  			x = float64(n)
	 349  			n <<= 2
	 350  		}
	 351  	}
	 352  	switch n {
	 353  	case 0:
	 354  		x = Sin(Pi * x)
	 355  	case 1, 2:
	 356  		x = Cos(Pi * (0.5 - x))
	 357  	case 3, 4:
	 358  		x = Sin(Pi * (1 - x))
	 359  	case 5, 6:
	 360  		x = -Cos(Pi * (x - 1.5))
	 361  	default:
	 362  		x = Sin(Pi * (x - 2))
	 363  	}
	 364  	return -x
	 365  }
	 366  

View as plain text