time-to-botec

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

higher_newton.js (3087B)

      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 'use strict';
     20 
     21 // MODULES //
     22 
     23 var logger = require( 'debug' );
     24 var gammainc = require( './../../../../base/special/gammainc' );
     25 var abs = require( './../../../../base/special/abs' );
     26 var exp = require( './../../../../base/special/exp' );
     27 var ln = require( './../../../../base/special/ln' );
     28 var MAX_FLOAT32 = require( '@stdlib/constants/float32/max' );
     29 
     30 
     31 // VARIABLES //
     32 
     33 var debug = logger( 'gammaincinv:higher_newton' );
     34 
     35 
     36 // MAIN //
     37 
     38 /**
     39 * Implementation of the high order Newton-like method.
     40 *
     41 * @private
     42 * @param {number} x0 - initial value
     43 * @param {number} a - scale parameter
     44 * @param {number} m - indicator
     45 * @param {Probability} p - probability value
     46 * @param {Probability} q - probability value
     47 * @param {number} lgama - logarithm of scale parameter
     48 * @param {number} invfp - one over `fp`
     49 * @param {boolean} pcase - boolean indicating whether p < 0.5
     50 * @returns {number} function value of the inverse
     51 */
     52 function higherNewton( x0, a, m, p, q, lgama, invfp, pcase ) {
     53 	var dlnr;
     54 	var xini;
     55 	var ck0;
     56 	var ck1;
     57 	var ck2;
     58 	var a2;
     59 	var x2;
     60 	var px;
     61 	var qx;
     62 	var xr;
     63 	var t;
     64 	var n;
     65 	var r;
     66 	var x;
     67 
     68 	x = x0;
     69 	t = 1;
     70 	n = 1;
     71 	a2 = a * a;
     72 	xini = x0;
     73 	do {
     74 		x = x0;
     75 		x2 = x * x;
     76 		if ( m === 0 ) {
     77 			dlnr = ( ( 1.0-a ) * ln( x ) ) + x + lgama;
     78 			if ( dlnr > ln( MAX_FLOAT32 ) ) {
     79 				debug( 'Warning: overflow problems in one or more steps of the computation. The initial approximation to the root is returned.' );
     80 				return xini;
     81 			}
     82 			r = exp( dlnr );
     83 		} else {
     84 			r = -invfp * x;
     85 		}
     86 		if ( pcase ) {
     87 			// Call: gammainc( x, s[, regularized = true ][, upper = false ] )
     88 			px = gammainc( x, a, true, false );
     89 			ck0 = -r * ( px - p );
     90 		} else {
     91 			// Call: gammainc( x, s[, regularized = true ][, upper = true ] )
     92 			qx = gammainc( x, a, true, true );
     93 			ck0 = r * ( qx - q );
     94 		}
     95 		r = ck0;
     96 		if ( ( p > 1e-120 ) || ( n > 1 ) ) {
     97 			ck1 = 0.5 * ( x - a + 1.0 ) / x;
     98 			ck2 = ( (2*x2) - (4*x*a) + (4*x) + (2*a2) - (3*a) + 1 ) / x2;
     99 			ck2 /= 6.0;
    100 			x0 = x + ( r * ( 1.0 + ( r * ( ck1 + (r*ck2) ) ) ) );
    101 		} else {
    102 			x0 = x + r;
    103 		}
    104 		t = abs( ( x/x0 ) - 1.0 );
    105 		n += 1;
    106 		x = x0;
    107 		if ( x < 0 ) {
    108 			x = xini;
    109 			n = 100;
    110 		}
    111 	} while ( ( ( t > 2e-14 ) && ( n < 35 ) ) );
    112 	if ( ( t > 2e-14 ) || ( n > 99 ) ) {
    113 		debug( 'Warning: the number of iterations in the Newton method reached the upper limit N=35. The last value obtained for the root is given as output.' );
    114 	}
    115 	xr = x || 0;
    116 	return xr;
    117 }
    118 
    119 
    120 // EXPORTS //
    121 
    122 module.exports = higherNewton;