time-to-botec

gammaln.js (10225B)

---

 /**
 * @license Apache-2.0
 *
 * Copyright (c) 2018 The Stdlib Authors.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *    http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 *
 * ## Notice
 *
 * The following copyright, license, and long comment were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/9.3.0/lib/msun/src/e_lgamma_r.c}. The implementation follows the original, but has been modified for JavaScript.
 *
 * ```text
 * 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.
 * ```
 */
 
 'use strict';
 
 // MODULES //
 
 var isnan = require( './../../../../base/assert/is-nan' );
 var isInfinite = require( './../../../../base/assert/is-infinite' );
 var abs = require( './../../../../base/special/abs' );
 var ln = require( './../../../../base/special/ln' );
 var trunc = require( './../../../../base/special/trunc' );
 var sinpi = require( './../../../../base/special/sinpi' );
 var PI = require( '@stdlib/constants/float64/pi' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var polyvalA1 = require( './polyval_a1.js' );
 var polyvalA2 = require( './polyval_a2.js' );
 var polyvalR = require( './polyval_r.js' );
 var polyvalS = require( './polyval_s.js' );
 var polyvalT1 = require( './polyval_t1.js' );
 var polyvalT2 = require( './polyval_t2.js' );
 var polyvalT3 = require( './polyval_t3.js' );
 var polyvalU = require( './polyval_u.js' );
 var polyvalV = require( './polyval_v.js' );
 var polyvalW = require( './polyval_w.js' );
 
 
 // VARIABLES //
 
 var A1C = 7.72156649015328655494e-02; // 0x3FB3C467E37DB0C8
 var A2C = 3.22467033424113591611e-01; // 0x3FD4A34CC4A60FAD
 var RC = 1.0;
 var SC = -7.72156649015328655494e-02; // 0xBFB3C467E37DB0C8
 var T1C = 4.83836122723810047042e-01; // 0x3FDEF72BC8EE38A2
 var T2C = -1.47587722994593911752e-01; // 0xBFC2E4278DC6C509
 var T3C = 6.46249402391333854778e-02; // 0x3FB08B4294D5419B
 var UC = -7.72156649015328655494e-02; // 0xBFB3C467E37DB0C8
 var VC = 1.0;
 var WC = 4.18938533204672725052e-01; // 0x3FDACFE390C97D69
 var YMIN = 1.461632144968362245;
 var TWO52 = 4503599627370496; // 2**52
 var TWO58 = 288230376151711744; // 2**58
 var TINY = 8.470329472543003e-22;
 var TC = 1.46163214496836224576e+00; // 0x3FF762D86356BE3F
 var TF = -1.21486290535849611461e-01; // 0xBFBF19B9BCC38A42
 var TT = -3.63867699703950536541e-18; // 0xBC50C7CAA48A971F => TT = -(tail of TF)
 
 
 // MAIN //
 
 /**
 * Evaluates the natural logarithm of the gamma function.
 *
 * ## Method
 *
 * 1.  Argument reduction for \\(0 < x \leq 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
 *
 *     ```tex
 *     \operatorname{lgamma}(1+s) = \ln(s) + \operatorname{lgamma}(s)
 *     ```
 *
 *     For example,
 *
 *     ```tex
 *     \begin{align*}
 *     \operatorname{lgamma}(7.3) &= \ln(6.3) + \operatorname{lgamma}(6.3) \\
 *     &= \ln(6.3 \cdot 5.3) + \operatorname{lgamma}(5.3) \\
 *     &= \ln(6.3 \cdot 5.3 \cdot 4.3 \cdot 3.3 \cdot2.3) + \operatorname{lgamma}(2.3)
 *     \end{align*}
 *     ```
 *
 * 2.  Compute a polynomial approximation of \\(\mathrm{lgamma}\\) around its minimum (\\(\mathrm{ymin} = 1.461632144968362245\\)) to maintain monotonicity. On the interval \\(\[\mathrm{ymin} - 0.23, \mathrm{ymin} + 0.27]\\) (i.e., \\(\[1.23164,1.73163]\\)), we let \\(z = x - \mathrm{ymin}\\) and use
 *
 *     ```tex
 *     \operatorname{lgamma}(x) = -1.214862905358496078218 + z^2 \cdot \operatorname{poly}(z)
 *     ```
 *
 *     where \\(\operatorname{poly}(z)\\) is a \\(14\\) degree polynomial.
 *
 * 3.  Compute a rational approximation in the primary interval \\(\[2,3]\\). Let \\( s = x - 2.0 \\). We can thus use the approximation
 *
 *     ```tex
 *     \operatorname{lgamma}(x) = \frac{s}{2} + s\frac{\operatorname{P}(s)}{\operatorname{Q}(s)}
 *     ```
 *
 *     with accuracy
 *
 *     ```tex
 *     \biggl|\frac{\mathrm{P}}{\mathrm{Q}} - \biggr(\operatorname{lgamma}(x)-\frac{s}{2}\biggl)\biggl| < 2^{-61.71}
 *     ```
 *
 *     The algorithms are based on the observation
 *
 *     ```tex
 *     \operatorname{lgamma}(2+s) = s(1 - \gamma) + \frac{\zeta(2) - 1}{2} s^2 - \frac{\zeta(3) - 1}{3} s^3 + \ldots
 *     ```
 *
 *     where \\(\zeta\\) is the zeta function and \\(\gamma = 0.5772156649...\\) is the Euler-Mascheroni constant, which is very close to \\(0.5\\).
 *
 * 4.  For \\(x \geq 8\\),
 *
 *     ```tex
 *     \operatorname{lgamma}(x) \approx \biggl(x-\frac{1}{2}\biggr) \ln(x) - x + \frac{\ln(2\pi)}{2} + \frac{1}{12x} - \frac{1}{360x^3} + \ldots
 *     ```
 *
 *     which can be expressed
 *
 *     ```tex
 *     \operatorname{lgamma}(x) \approx \biggl(x-\frac{1}{2}\biggr)(\ln(x)-1)-\frac{\ln(2\pi)-1}{2} + \ldots
 *     ```
 *
 *     Let \\(z = \frac{1}{x}\\). We can then use the approximation
 *
 *     ```tex
 *     f(z) = \operatorname{lgamma}(x) - \biggl(x-\frac{1}{2}\biggr)(\ln(x)-1)
 *     ```
 *
 *     by
 *
 *     ```tex
 *     w = w_0 + w_1 z + w_2 z^3 + w_3 z^5 + \ldots + w_6 z^{11}
 *     ```
 *
 *     where
 *
 *     ```tex
 *     |w - f(z)| < 2^{-58.74}
 *     ```
 *
 * 5.  For negative \\(x\\), since
 *
 *     ```tex
 *     -x \Gamma(-x) \Gamma(x) = \frac{\pi}{\sin(\pi x)}
 *     ```
 *
 *     where \\(\Gamma\\) is the gamma function, we have
 *
 *     ```tex
 *     \Gamma(x) = \frac{\pi}{\sin(\pi x)(-x)\Gamma(-x)}
 *     ```
 *
 *     Since \\(\Gamma(-x)\\) is positive,
 *
 *     ```tex
 *     \operatorname{sign}(\Gamma(x)) = \operatorname{sign}(\sin(\pi x))
 *     ```
 *
 *     for \\(x < 0\\). Hence, for \\(x < 0\\),
 *
 *     ```tex
 *     \mathrm{signgam} = \operatorname{sign}(\sin(\pi x))
 *     ```
 *
 *     and
 *
 *     ```tex
 *     \begin{align*}
 *     \operatorname{lgamma}(x) &= \ln(|\Gamma(x)|) \\
 *     &= \ln\biggl(\frac{\pi}{|x \sin(\pi x)|}\biggr) - \operatorname{lgamma}(-x)
 *     \end{align*}
 *     ```
 *
 *     <!-- <note> -->
 *
 *     Note that one should avoid computing \\(\pi (-x)\\) directly in the computation of \\(\sin(\pi (-x))\\).
 *
 *     <!-- </note> -->
 *
 *
 * ## Special Cases
 *
 * ```tex
 * \begin{align*}
 * \operatorname{lgamma}(2+s) &\approx s (1-\gamma) & \mathrm{for\ tiny\ s} \\
 * \operatorname{lgamma}(x) &\approx -\ln(x) & \mathrm{for\ tiny\ x} \\
 * \operatorname{lgamma}(1) &= 0 & \\
 * \operatorname{lgamma}(2) &= 0 & \\
 * \operatorname{lgamma}(0) &= \infty & \\
 * \operatorname{lgamma}(\infty) &= \infty & \\
 * \operatorname{lgamma}(-\mathrm{integer}) &= \pm \infty
 * \end{align*}
 * ```
 *
 *
 * @param {number} x - input value
 * @returns {number} function value
 *
 * @example
 * var v = gammaln( 1.0 );
 * // returns 0.0
 *
 * @example
 * var v = gammaln( 2.0 );
 * // returns 0.0
 *
 * @example
 * var v = gammaln( 4.0 );
 * // returns ~1.792
 *
 * @example
 * var v = gammaln( -0.5 );
 * // returns ~1.266
 *
 * @example
 * var v = gammaln( 0.5 );
 * // returns ~0.572
 *
 * @example
 * var v = gammaln( 0.0 );
 * // returns Infinity
 *
 * @example
 * var v = gammaln( NaN );
 * // returns NaN
 */
 function gammaln( x ) {
 	var isNegative;
 	var nadj;
 	var flg;
 	var p3;
 	var p2;
 	var p1;
 	var p;
 	var q;
 	var t;
 	var w;
 	var y;
 	var z;
 	var r;
 
 	// Special cases: NaN, +-infinity
 	if ( isnan( x ) || isInfinite( x ) ) {
 		return x;
 	}
 	// Special case: 0
 	if ( x === 0.0 ) {
 		return PINF;
 	}
 	if ( x < 0.0 ) {
 		isNegative = true;
 		x = -x;
 	} else {
 		isNegative = false;
 	}
 	// If |x| < 2**-70, return -ln(|x|)
 	if ( x < TINY ) {
 		return -ln( x );
 	}
 	if ( isNegative ) {
 		// If |x| >= 2**52, must be -integer
 		if ( x >= TWO52 ) {
 			return PINF;
 		}
 		t = sinpi( x );
 		if ( t === 0.0 ) {
 			return PINF;
 		}
 		nadj = ln( PI / abs( t*x ) );
 	}
 	// If x equals 1 or 2, return 0
 	if ( x === 1.0 || x === 2.0 ) {
 		return 0.0;
 	}
 	// If x < 2, use lgamma(x) = lgamma(x+1) - log(x)
 	if ( x < 2.0 ) {
 		if ( x <= 0.9 ) {
 			r = -ln( x );
 
 			// 0.7316 <= x <=  0.9
 			if ( x >= ( YMIN - 1.0 + 0.27 ) ) {
 				y = 1.0 - x;
 				flg = 0;
 			}
 			// 0.2316 <= x < 0.7316
 			else if ( x >= (YMIN - 1.0 - 0.27) ) {
 				y = x - (TC - 1.0);
 				flg = 1;
 			}
 			// 0 < x < 0.2316
 			else {
 				y = x;
 				flg = 2;
 			}
 		} else {
 			r = 0.0;
 
 			// 1.7316 <= x < 2
 			if ( x >= (YMIN + 0.27) ) {
 				y = 2.0 - x;
 				flg = 0;
 			}
 			// 1.2316 <= x < 1.7316
 			else if ( x >= (YMIN - 0.27) ) {
 				y = x - TC;
 				flg = 1;
 			}
 			// 0.9 < x < 1.2316
 			else {
 				y = x - 1.0;
 				flg = 2;
 			}
 		}
 		switch ( flg ) { // eslint-disable-line default-case
 		case 0:
 			z = y * y;
 			p1 = A1C + (z*polyvalA1( z ));
 			p2 = z * (A2C + (z*polyvalA2( z )));
 			p = (y*p1) + p2;
 			r += ( p - (0.5*y) );
 			break;
 		case 1:
 			z = y * y;
 			w = z * y;
 			p1 = T1C + (w*polyvalT1( w ));
 			p2 = T2C + (w*polyvalT2( w ));
 			p3 = T3C + (w*polyvalT3( w ));
 			p = (z*p1) - (TT - (w*(p2+(y*p3))));
 			r += ( TF + p );
 			break;
 		case 2:
 			p1 = y * (UC + (y*polyvalU( y )));
 			p2 = VC + (y*polyvalV( y ));
 			r += (-0.5*y) + (p1/p2);
 			break;
 		}
 	}
 	// 2 <= x < 8
 	else if ( x < 8.0 ) {
 		flg = trunc( x );
 		y = x - flg;
 		p = y * (SC + (y*polyvalS( y )));
 		q = RC + (y*polyvalR( y ));
 		r = (0.5*y) + (p/q);
 		z = 1.0; // gammaln(1+s) = ln(s) + gammaln(s)
 		switch ( flg ) { // eslint-disable-line default-case
 		case 7:
 			z *= y + 6.0;
 
 			/* falls through */
 		case 6:
 			z *= y + 5.0;
 
 			/* falls through */
 		case 5:
 			z *= y + 4.0;
 
 			/* falls through */
 		case 4:
 			z *= y + 3.0;
 
 			/* falls through */
 		case 3:
 			z *= y + 2.0;
 			r += ln( z );
 		}
 	}
 	// 8 <= x < 2**58
 	else if ( x < TWO58 ) {
 		t = ln( x );
 		z = 1.0 / x;
 		y = z * z;
 		w = WC + (z*polyvalW( y ));
 		r = ((x-0.5)*(t-1.0)) + w;
 	}
 	// 2**58 <= x <= Inf
 	else {
 		r = x * ( ln(x)-1.0 );
 	}
 	if ( isNegative ) {
 		r = nadj - r;
 	}
 	return r;
 }
 
 
 // EXPORTS //
 
 module.exports = gammaln;