time-to-botec

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

beta.js (3568B)

      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]{http://www.boost.org/doc/libs/1_64_0/boost/math/special_functions/beta.hpp}. The implementation has been modified for JavaScript.
     22 *
     23 * ```text
     24 * (C) Copyright John Maddock 2006.
     25 *
     26 * Use, modification and distribution are subject to the
     27 * Boost Software License, Version 1.0. (See accompanying file
     28 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
     29 * ```
     30 */
     31 
     32 'use strict';
     33 
     34 // MODULES //
     35 
     36 var isnan = require( './../../../../base/assert/is-nan' );
     37 var log1p = require( './../../../../base/special/log1p' );
     38 var sqrt = require( './../../../../base/special/sqrt' );
     39 var abs = require( './../../../../base/special/abs' );
     40 var exp = require( './../../../../base/special/exp' );
     41 var pow = require( './../../../../base/special/pow' );
     42 var E = require( '@stdlib/constants/float64/e' );
     43 var EPSILON = require( '@stdlib/constants/float64/eps' );
     44 var lanczosSumExpGScaled = require( './lanczos_sum_expg_scaled.js' ); // Lanczos approximation scaled by exp(G)
     45 
     46 
     47 // VARIABLES //
     48 
     49 var G = 10.90051099999999983936049829935654997826;
     50 
     51 
     52 // MAIN //
     53 
     54 /**
     55 * Evaluate the beta function.
     56 *
     57 * @param {NonNegativeNumber} a - input value
     58 * @param {NonNegativeNumber} b - input value
     59 * @returns {number} evaluated beta function
     60 *
     61 * @example
     62 * var v = beta( 0.0, 0.5 );
     63 * // returns Infinity
     64 *
     65 * @example
     66 * var v = beta( 1.0, 1.0 );
     67 * // returns 1.0
     68 *
     69 * @example
     70 * var v = beta( -1.0, 2.0 );
     71 * // returns NaN
     72 *
     73 * @example
     74 * var v = beta( 5.0, 0.2 );
     75 * // returns ~3.382
     76 *
     77 * @example
     78 * var v = beta( 4.0, 1.0 );
     79 * // returns 0.25
     80 *
     81 * @example
     82 * var v = beta( NaN, 2.0 );
     83 * // returns NaN
     84 */
     85 function beta( a, b ) {
     86 	var ambh;
     87 	var agh;
     88 	var bgh;
     89 	var cgh;
     90 	var res;
     91 	var tmp;
     92 	var c;
     93 
     94 	if ( isnan( a ) || isnan( b ) ) {
     95 		return NaN;
     96 	}
     97 	if ( a < 0.0 || b < 0.0 ) {
     98 		return NaN;
     99 	}
    100 	if ( b === 1.0 ) {
    101 		return 1.0 / a;
    102 	}
    103 	if ( a === 1.0 ) {
    104 		return 1.0 / b;
    105 	}
    106 	c = a + b;
    107 	if ( c < EPSILON ) {
    108 		res = c / a;
    109 		res /= b;
    110 		return res;
    111 	}
    112 
    113 	// Special cases:
    114 	if ( c === a && b < EPSILON ) {
    115 		return 1.0 / b;
    116 	}
    117 	if ( c === b && a < EPSILON ) {
    118 		return 1.0 / a;
    119 	}
    120 
    121 	if ( a < b ) {
    122 		// Swap `a` and `b`:
    123 		tmp = b;
    124 		b = a;
    125 		a = tmp;
    126 	}
    127 
    128 	// Lanczos calculation:
    129 	agh = a + G - 0.5;
    130 	bgh = b + G - 0.5;
    131 	cgh = c + G - 0.5;
    132 	res = lanczosSumExpGScaled( a ) * ( lanczosSumExpGScaled( b )/lanczosSumExpGScaled( c ) ); // eslint-disable-line max-len
    133 	ambh = a - 0.5 - b;
    134 	if ( ( abs( b*ambh ) < ( cgh*100.0 ) ) && a > 100.0 ) {
    135 		// Special case where the base of the power term is close to 1; compute `(1+x)^y` instead:
    136 		res *= exp( ambh * log1p( -b/cgh ) );
    137 	} else {
    138 		res *= pow( agh/cgh, ambh );
    139 	}
    140 	if ( cgh > 1.0e10 ) {
    141 		// This avoids possible overflow, but appears to be marginally less accurate:
    142 		res *= pow( (agh/cgh)*(bgh/cgh), b );
    143 	} else {
    144 		res *= pow( (agh*bgh)/(cgh*cgh), b );
    145 	}
    146 	res *= sqrt( E/bgh);
    147 	return res;
    148 }
    149 
    150 
    151 // EXPORTS //
    152 
    153 module.exports = beta;