time-to-botec

Benchmark sampling in different programming languages
Log | Files | Refs | README

gamma.js (4452B)

      1 /**
      2 * @license Apache-2.0
      3 *
      4 * Copyright (c) 2018 The Stdlib Authors.
      5 *
      6 * Licensed under the Apache License, Version 2.0 (the "License");
      7 * you may not use this file except in compliance with the License.
      8 * You may obtain a copy of the License at
      9 *
     10 *    http://www.apache.org/licenses/LICENSE-2.0
     11 *
     12 * Unless required by applicable law or agreed to in writing, software
     13 * distributed under the License is distributed on an "AS IS" BASIS,
     14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     15 * See the License for the specific language governing permissions and
     16 * limitations under the License.
     17 *
     18 *
     19 * ## Notice
     20 *
     21 * The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_64_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
     22 *
     23 * ```text
     24 * Copyright John Maddock 2006-7, 2013-14.
     25 * Copyright Paul A. Bristow 2007, 2013-14.
     26 * Copyright Nikhar Agrawal 2013-14.
     27 * Copyright Christopher Kormanyos 2013-14.
     28 *
     29 * Use, modification and distribution are subject to the
     30 * Boost Software License, Version 1.0. (See accompanying file
     31 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
     32 * ```
     33 */
     34 
     35 'use strict';
     36 
     37 // MODULES //
     38 
     39 var gammaLanczosSum = require( './../../../../../base/special/gamma-lanczos-sum' );
     40 var isNegativeZero = require( './../../../../../base/assert/is-negative-zero' );
     41 var isInteger = require( './../../../../../base/assert/is-integer' );
     42 var isnan = require( './../../../../../base/assert/is-nan' );
     43 var signum = require( './../../../../../base/special/signum' );
     44 var floor = require( './../../../../../base/special/floor' );
     45 var abs = require( './../../../../../base/special/abs' );
     46 var exp = require( './../../../../../base/special/exp' );
     47 var pow = require( './../../../../../base/special/pow' );
     48 var ln = require( './../../../../../base/special/ln' );
     49 var SQRT_EPSILON = require( '@stdlib/constants/float64/sqrt-eps' );
     50 var LANCZOS_G = require( '@stdlib/constants/float64/gamma-lanczos-g' );
     51 var EULERGAMMA = require( '@stdlib/constants/float64/eulergamma' );
     52 var MAX_VALUE = require( '@stdlib/constants/float64/max' );
     53 var MAX_LN = require( '@stdlib/constants/float64/max-ln' );
     54 var PINF = require( '@stdlib/constants/float64/pinf' );
     55 var NINF = require( '@stdlib/constants/float64/ninf' );
     56 var PI = require( '@stdlib/constants/float64/pi' );
     57 var sinpx = require( './sinpx.js' );
     58 var FACTORIALS = require( './factorials.json' );
     59 
     60 
     61 // VARIABLES //
     62 
     63 var MAX_FACTORIAL = 170;
     64 var HALF = 0.5;
     65 
     66 
     67 // MAIN //
     68 
     69 /**
     70 * Evaluates the gamma function.
     71 *
     72 * @param {number} x - input value
     73 * @returns {number} function value
     74 *
     75 * @example
     76 * var v = gamma( 4.0 );
     77 * // returns 6.0
     78 *
     79 * @example
     80 * var v = gamma( -1.5 );
     81 * // returns ~2.363
     82 *
     83 * @example
     84 * var v = gamma( -0.5 );
     85 * // returns ~-3.545
     86 *
     87 * @example
     88 * var v = gamma( 0.5 );
     89 * // returns ~1.772
     90 *
     91 * @example
     92 * var v = gamma( 0.0 );
     93 * // returns Infinity
     94 *
     95 * @example
     96 * var v = gamma( -0.0 );
     97 * // returns -Infinity
     98 *
     99 * @example
    100 * var v = gamma( NaN );
    101 * // returns NaN
    102 */
    103 function gamma( x ) {
    104 	var result;
    105 	var lzgh;
    106 	var zgh;
    107 	var hp;
    108 
    109 	if ( x === NINF || isnan( x ) ) {
    110 		return NaN;
    111 	}
    112 	if ( x === 0.0 ) {
    113 		if ( isNegativeZero( x ) ) {
    114 			return NINF;
    115 		}
    116 		return PINF;
    117 	}
    118 	result = 1.0;
    119 	if ( x < 0.0 ) {
    120 		if ( isInteger( x ) ) {
    121 			return NaN;
    122 		}
    123 		if ( x <= -20.0 ) {
    124 			result = gamma( -x ) * sinpx( x );
    125 			if ( abs(result) < 1.0 && MAX_VALUE * abs(result) < PI ) {
    126 				return ( signum( result ) === 1 ) ? NINF : PINF;
    127 			}
    128 			result = -PI / result;
    129 			return result;
    130 		}
    131 		// Shift x to > 1:
    132 		while ( x < 0.0 ) {
    133 			result /= x;
    134 			x += 1.0;
    135 		}
    136 	}
    137 	if ( floor( x ) === x && x < MAX_FACTORIAL ) {
    138 		result *= FACTORIALS[ floor( x ) - 1 ];
    139 	}
    140 	else if ( x < SQRT_EPSILON ) {
    141 		if ( x < 1.0 / MAX_VALUE ) {
    142 			result = PINF;
    143 		}
    144 		result *= ( 1.0 / x ) - EULERGAMMA;
    145 	}
    146 	else {
    147 		result *= gammaLanczosSum( x );
    148 		zgh = ( x + LANCZOS_G - HALF );
    149 		lzgh = ln( zgh );
    150 		if ( x * lzgh > MAX_LN ) {
    151 			// We're going to overflow unless this is done with care:
    152 			if ( lzgh * x / 2.0 > MAX_LN ) {
    153 				return ( signum( result ) === 1 ) ? PINF : NINF;
    154 			}
    155 			hp = pow( zgh, ( x / 2.0 ) - 0.25 );
    156 			result *= hp / exp( zgh );
    157 			if ( MAX_VALUE / hp < result ) {
    158 				return ( signum( result ) === 1 ) ? PINF : NINF;
    159 			}
    160 			result *= hp;
    161 		} else {
    162 			result *= pow( zgh, x - HALF ) / exp( zgh );
    163 		}
    164 	}
    165 	return result;
    166 }
    167 
    168 
    169 // EXPORTS //
    170 
    171 module.exports = gamma;