time-to-botec

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

lgamma_small_imp.js (4422B)

      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/detail/lgamma_small.hpp}. The implementation has been modified for JavaScript.
     22 *
     23 * ```text
     24 * (C) Copyright John Maddock 2006-7, 2013-14.
     25 * (C) Copyright Paul A. Bristow 2007, 2013-14.
     26 * (C) Copyright Nikhar Agrawal 2013-14.
     27 * (C) 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 ln = require( './../../../../base/special/ln' );
     40 var EPS = require( '@stdlib/constants/float64/eps' );
     41 var rateval1 = require( './rational_p1q1.js' );
     42 var rateval2 = require( './rational_p2q2.js' );
     43 var rateval3 = require( './rational_p3q3.js' );
     44 
     45 
     46 // VARIABLES //
     47 
     48 var Y1 = 0.158963680267333984375;
     49 var Y2 = 0.52815341949462890625;
     50 var Y3 = 0.452017307281494140625;
     51 
     52 
     53 // MAIN //
     54 
     55 /**
     56 * Evaluates the natural logarithm of the gamma function for small arguments.
     57 *
     58 * ## Method
     59 *
     60 * 1.  For \\( z > 2 \\), begin by performing argument reduction until \\( z \\) is in \\(\[2,3)\\). Use the following form:
     61 *
     62 *     ```tex
     63 *     \operatorname{gammaln}(z) = (z-2)(z+1)(Y + R(z-2))
     64 *     ```
     65 *
     66 *     where \\( R(z-2) \\) is a rational approximation optimized for low absolute error. As long as the absolute error is small compared to the constant \\( Y \\), then any rounding error in the computation will get wiped out.
     67 *
     68 * 2.  If \\( z < 1 \\), use recurrence to shift to \\( z \\) in the interval \\(\[1,2\]\\). Then, use one of two approximations: one for \\( z \\) in \\(\[1,1.5\]\\) and one for \\( z \\) in \\(\[1.5,2\]\\):
     69 *
     70 *     -   For \(( z \\) in \\(\[1,1.5\]\\), use
     71 *
     72 *         ```tex
     73 *         \operatorname{gammaln}(z) = (z-1)(z-2)(Y + R(z-1))
     74 *         ```
     75 *
     76 *         where \\( R(z-1) \\) is a rational approximation optimized for low absolute error. As long as the absolute error is small compared to the constant \\( Y \\), then any rounding error in the computation will get wiped out.
     77 *
     78 *     -   For \\( z \\) in \\(\[1.5,2\]\\), use
     79 *
     80 *         ```tex
     81 *         \operatorname{gammaln}(z) = (2-z)(1-z)(Y + R(2-z))
     82 *         ```
     83 *
     84 *         where \\( R(2-z) \\) is a rational approximation optimized for low absolute error. As long as the absolute error is small compared to the constant \\( Y \\), then any rounding error in the computation will get wiped out.
     85 *
     86 *
     87 * ## Notes
     88 *
     89 * -   Relative error:
     90 *
     91 *     | function | peak         | maximum deviation |
     92 *     |:--------:|:------------:|:-----------------:|
     93 *     | R(Z-2)   | 4.231e-18    | 5.900e-24         |
     94 *     | R(Z-1)   | 1.230011e-17 | 3.139e-021        |
     95 *     | R(2-Z)   | 1.797565e-17 | 2.151e-021        |
     96 *
     97 *
     98 * @private
     99 * @param {number} z - input value
    100 * @param {number} zm1 - `z` minus one
    101 * @param {number} zm2 - `z` minus two
    102 * @returns {number} function value
    103 */
    104 function lgammaSmallImp( z, zm1, zm2 ) {
    105 	var prefix;
    106 	var result;
    107 	var r;
    108 	var R;
    109 
    110 	if ( z < EPS ) {
    111 		return -ln( z );
    112 	}
    113 	if ( zm1 === 0.0 || zm2 === 0.0 ) {
    114 		return 0.0;
    115 	}
    116 	result = 0.0;
    117 	if ( z > 2.0 ) {
    118 		if ( z >= 3.0 ) {
    119 			do {
    120 				z -= 1.0;
    121 				zm2 -= 1.0;
    122 				result += ln(z);
    123 			} while ( z >= 3.0 );
    124 			zm2 = z - 2.0;
    125 		}
    126 		r = zm2 * ( z+1.0 );
    127 		R = rateval1( zm2 );
    128 		result += ( r*Y1 ) + ( r*R );
    129 		return result;
    130 	}
    131 	if ( z < 1.0 ) {
    132 		result += -ln(z);
    133 		zm2 = zm1;
    134 		zm1 = z;
    135 		z += 1.0;
    136 	}
    137 	if ( z <= 1.5 ) {
    138 		r = rateval2( zm1 );
    139 		prefix = zm1 * zm2;
    140 		result += ( prefix*Y2 ) + ( prefix*r );
    141 		return result;
    142 	}
    143 	// Case: 1.5 < z <= 2
    144 	r = zm2 * zm1;
    145 	R = rateval3( -zm2 );
    146 	result += ( r*Y3 ) + ( r*R );
    147 	return result;
    148 }
    149 
    150 
    151 // EXPORTS //
    152 
    153 module.exports = lgammaSmallImp;