time-to-botec

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

regularised_gamma_prefix.js (4261B)

      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 lanczosSumExpGScaled = require( './../../../../base/special/gamma-lanczos-sum-expg-scaled' );
     40 var gammaln = require( './../../../../base/special/gammaln' );
     41 var gamma = require( './../../../../base/special/gamma' );
     42 var log1p = require( './../../../../base/special/log1p' );
     43 var sqrt = require( './../../../../base/special/sqrt' );
     44 var abs = require( './../../../../base/special/abs' );
     45 var exp = require( './../../../../base/special/exp' );
     46 var pow = require( './../../../../base/special/pow' );
     47 var max = require( './../../../../base/special/max' );
     48 var min = require( './../../../../base/special/min' );
     49 var ln = require( './../../../../base/special/ln' );
     50 var MAX_LN = require( '@stdlib/constants/float64/max-ln' );
     51 var MIN_LN = require( '@stdlib/constants/float64/min-ln' );
     52 var G = require( '@stdlib/constants/float64/gamma-lanczos-g' );
     53 var E = require( '@stdlib/constants/float64/e' );
     54 
     55 
     56 // MAIN //
     57 
     58 /**
     59 * Computes `(z^a)*(e^-z) / gamma(a)`.
     60 *
     61 * @private
     62 * @param {number} a - input value
     63 * @param {number} z - input value
     64 * @returns {number} function value
     65 */
     66 function regularisedGammaPrefix( a, z ) {
     67 	var prefix;
     68 	var amza;
     69 	var agh;
     70 	var alz;
     71 	var amz;
     72 	var sq;
     73 	var d;
     74 
     75 	agh = a + G - 0.5;
     76 	d = ( (z - a) - G + 0.5 ) / agh;
     77 	if ( a < 1.0 ) {
     78 		// Treat a < 1 as a special case because our Lanczos approximations are optimized against the factorials with a > 1, and for high precision types very small values of `a` can give rather erroneous results for gamma:
     79 		if ( z <= MIN_LN ) {
     80 			// Use logs, so should be free of cancellation errors:
     81 			return exp( ( a * ln(z) ) - z - gammaln( a ) );
     82 		}
     83 		// No danger of overflow as gamma(a) < 1/a for small a, so direct calculation:
     84 		return pow( z, a ) * exp( -z ) / gamma( a );
     85 	}
     86 	if ( abs(d*d*a) <= 100.0 && a > 150.0 ) {
     87 		// Special case for large a and a ~ z:
     88 		prefix = ( a * ( log1p( d ) - d ) ) + ( z * ( 0.5-G ) / agh );
     89 		prefix = exp( prefix );
     90 	}
     91 	else {
     92 		// General case. Direct computation is most accurate, but use various fallbacks for different parts of the problem domain:
     93 		alz = a * ln(z / agh);
     94 		amz = a - z;
     95 		if (
     96 			min(alz, amz) <= MIN_LN ||
     97 			max(alz, amz) >= MAX_LN
     98 		) {
     99 			amza = amz / a;
    100 			if (
    101 				min(alz, amz)/2.0 > MIN_LN &&
    102 				max(alz, amz)/2.0 < MAX_LN
    103 			) {
    104 				// Compute square root of the result and then square it:
    105 				sq = pow( z / agh, a / 2.0 ) * exp( amz / 2.0 );
    106 				prefix = sq * sq;
    107 			}
    108 			else if (
    109 				min(alz, amz)/4.0 > MIN_LN &&
    110 				max(alz, amz)/4.0 < MAX_LN &&
    111 				z > a
    112 			) {
    113 				// Compute the 4th root of the result then square it twice:
    114 				sq = pow( z / agh, a / 4.0 ) * exp( amz / 4.0 );
    115 				prefix = sq * sq;
    116 				prefix *= prefix;
    117 			}
    118 			else if (
    119 				amza > MIN_LN &&
    120 				amza < MAX_LN
    121 			) {
    122 				prefix = pow( (z * exp(amza)) / agh, a );
    123 			}
    124 			else {
    125 				prefix = exp( alz + amz );
    126 			}
    127 		}
    128 		else
    129 		{
    130 			prefix = pow( z / agh, a ) * exp( amz );
    131 		}
    132 	}
    133 	prefix *= sqrt( agh / E ) / lanczosSumExpGScaled( a );
    134 	return prefix;
    135 }
    136 
    137 
    138 // EXPORTS //
    139 
    140 module.exports = regularisedGammaPrefix;