time-to-botec

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

gamma.js (4531B)

      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, long comment, copyright, license, and constants are from [Cephes]{@link http://www.netlib.org/cephes}. The implementation follows the original, but has been modified for JavaScript.
     22 *
     23 * ```text
     24 * Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
     25 *
     26 * 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.
     27 *
     28 * Stephen L. Moshier
     29 * moshier@na-net.ornl.gov
     30 * ```
     31 */
     32 
     33 'use strict';
     34 
     35 // MODULES //
     36 
     37 var isnan = require( './../../../../base/assert/is-nan' );
     38 var isInteger = require( './../../../../base/assert/is-integer' );
     39 var isNegativeZero = require( './../../../../base/assert/is-negative-zero' );
     40 var abs = require( './../../../../base/special/abs' );
     41 var floor = require( './../../../../base/special/floor' );
     42 var sin = require( './../../../../base/special/sin' );
     43 var PINF = require( '@stdlib/constants/float64/pinf' );
     44 var NINF = require( '@stdlib/constants/float64/ninf' );
     45 var PI = require( '@stdlib/constants/float64/pi' );
     46 var stirlingApprox = require( './stirling_approximation.js' );
     47 var smallApprox = require( './small_approximation.js' );
     48 var rateval = require( './rational_pq.js' );
     49 
     50 
     51 // MAIN //
     52 
     53 /**
     54 * Evaluates the gamma function.
     55 *
     56 * ## Method
     57 *
     58 * 1.  Arguments \\(|x| \leq 34\\) are reduced by recurrence and the function approximated by a rational function of degree \\(6/7\\) in the interval \\((2,3)\\).
     59 * 2.  Large negative arguments are made positive using a reflection formula.
     60 * 3.  Large arguments are handled by Stirling's formula.
     61 *
     62 *
     63 * ## Notes
     64 *
     65 * -   Relative error:
     66 *
     67 *     | arithmetic | domain    | # trials | peak    | rms     |
     68 *     |:----------:|:---------:|:--------:|:-------:|:-------:|
     69 *     | DEC        | -34,34    | 10000    | 1.3e-16 | 2.5e-17 |
     70 *     | IEEE       | -170,-33  | 20000    | 2.3e-15 | 3.3e-16 |
     71 *     | IEEE       | -33, 33   | 20000    | 9.4e-16 | 2.2e-16 |
     72 *     | IEEE       | 33, 171.6 | 20000    | 2.3e-15 | 3.2e-16 |
     73 *
     74 * -   Error for arguments outside the test range will be larger owing to error amplification by the exponential function.
     75 *
     76 *
     77 * @param {number} x - input value
     78 * @returns {number} function value
     79 *
     80 * @example
     81 * var v = gamma( 4.0 );
     82 * // returns 6.0
     83 *
     84 * @example
     85 * var v = gamma( -1.5 );
     86 * // returns ~2.363
     87 *
     88 * @example
     89 * var v = gamma( -0.5 );
     90 * // returns ~-3.545
     91 *
     92 * @example
     93 * var v = gamma( 0.5 );
     94 * // returns ~1.772
     95 *
     96 * @example
     97 * var v = gamma( 0.0 );
     98 * // returns Infinity
     99 *
    100 * @example
    101 * var v = gamma( -0.0 );
    102 * // returns -Infinity
    103 *
    104 * @example
    105 * var v = gamma( NaN );
    106 * // returns NaN
    107 */
    108 function gamma( x ) {
    109 	var sign;
    110 	var q;
    111 	var p;
    112 	var z;
    113 	if (
    114 		(isInteger( x ) && x < 0) ||
    115 		x === NINF ||
    116 		isnan( x )
    117 	) {
    118 		return NaN;
    119 	}
    120 	if ( x === 0.0 ) {
    121 		if ( isNegativeZero( x ) ) {
    122 			return NINF;
    123 		}
    124 		return PINF;
    125 	}
    126 	if ( x > 171.61447887182298 ) {
    127 		return PINF;
    128 	}
    129 	if ( x < -170.5674972726612 ) {
    130 		return 0.0;
    131 	}
    132 	q = abs( x );
    133 	if ( q > 33.0 ) {
    134 		if ( x >= 0.0 ) {
    135 			return stirlingApprox( x );
    136 		}
    137 		p = floor( q );
    138 
    139 		// Check whether `x` is even...
    140 		if ( (p&1) === 0 ) {
    141 			sign = -1.0;
    142 		} else {
    143 			sign = 1.0;
    144 		}
    145 		z = q - p;
    146 		if ( z > 0.5 ) {
    147 			p += 1.0;
    148 			z = q - p;
    149 		}
    150 		z = q * sin( PI * z );
    151 		return sign * PI / ( abs(z)*stirlingApprox(q) );
    152 	}
    153 	// Reduce `x`...
    154 	z = 1.0;
    155 	while ( x >= 3.0 ) {
    156 		x -= 1.0;
    157 		z *= x;
    158 	}
    159 	while ( x < 0.0 ) {
    160 		if ( x > -1.0e-9 ) {
    161 			return smallApprox( x, z );
    162 		}
    163 		z /= x;
    164 		x += 1.0;
    165 	}
    166 	while ( x < 2.0 ) {
    167 		if ( x < 1.0e-9 ) {
    168 			return smallApprox( x, z );
    169 		}
    170 		z /= x;
    171 		x += 1.0;
    172 	}
    173 	if ( x === 2.0 ) {
    174 		return z;
    175 	}
    176 	x -= 2.0;
    177 	return z * rateval( x );
    178 }
    179 
    180 
    181 // EXPORTS //
    182 
    183 module.exports = gamma;