time-to-botec

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

gamma_delta_ratio.js (3264B)

      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 abs = require( './../../../../base/special/abs' );
     40 var floor = require( './../../../../base/special/floor' );
     41 var gamma = require( './../../../../base/special/gamma' );
     42 var factorial = require( './../../../../base/special/factorial' );
     43 var gammaDeltaRatioLanczos = require( './gamma_delta_ratio_lanczos.js' );
     44 
     45 
     46 // VARIABLES //
     47 
     48 var MAX_FACTORIAL = 170; // TODO: consider moving to pkg
     49 
     50 
     51 // MAIN //
     52 
     53 /**
     54 * Computes the ratio of two gamma functions.
     55 *
     56 * ## Notes
     57 *
     58 * -   Specifically, the function evaluates
     59 *
     60 *     ```tex
     61 *     \frac{ \Gamma( z ) }{ \Gamma( z + \delta ) }
     62 *     ```
     63 *
     64 * @param {number} z - first gamma parameter
     65 * @param {number} delta - difference
     66 * @returns {number} gamma ratio
     67 *
     68 * @example
     69 * var y = gammaDeltaRatio( 2.0, 3.0 );
     70 * // returns ~0.042
     71 *
     72 * @example
     73 * var y = gammaDeltaRatio( 4.0, 0.5 );
     74 * // returns ~0.516
     75 *
     76 * @example
     77 * var y = gammaDeltaRatio( 100.0, 0.0 );
     78 * // returns 1.0
     79 */
     80 function gammaDeltaRatio( z, delta ) {
     81 	var result;
     82 	var idelta;
     83 	var iz;
     84 
     85 	if ( z <= 0.0 || z + delta <= 0.0 ) {
     86 		// This isn't very sophisticated, or accurate, but it does work:
     87 		return gamma( z ) / gamma( z + delta );
     88 	}
     89 	idelta = floor( delta );
     90 	if ( idelta === delta ) {
     91 		iz = floor( z );
     92 		if ( iz === z ) {
     93 			// As both `z` and `delta` are integers, see if we can use a table lookup:
     94 			if ( z <= MAX_FACTORIAL && ( z + delta <= MAX_FACTORIAL ) ) {
     95 				return factorial( iz - 1.0 ) / factorial( idelta + iz - 1.0 );
     96 			}
     97 		}
     98 		if ( abs(delta) < 20.0 ) {
     99 			// As `delta` is a small integer, we can use a finite product:
    100 			if ( delta === 0.0 ) {
    101 				return 1.0;
    102 			}
    103 			if ( delta < 0.0 ) {
    104 				z -= 1.0;
    105 				result = z;
    106 				delta += 1.0;
    107 				while ( delta !== 0.0 ) {
    108 					z -= 1.0;
    109 					result *= z;
    110 					delta += 1.0;
    111 				}
    112 				return result;
    113 			}
    114 			result = 1.0 / z;
    115 			delta -= 1.0;
    116 			while ( delta !== 0.0 ) {
    117 				z += 1.0;
    118 				result /= z;
    119 				delta -= 1.0;
    120 			}
    121 			return result;
    122 		}
    123 	}
    124 	return gammaDeltaRatioLanczos( z, delta );
    125 }
    126 
    127 
    128 // EXPORTS //
    129 
    130 module.exports = gammaDeltaRatio;