time-to-botec

temme2.js (7687B)

---

 /**
 * @license Apache-2.0
 *
 * Copyright (c) 2018 The Stdlib Authors.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *    http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 *
 * ## Notice
 *
 * 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/ibeta_inverse.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * Copyright John Maddock 2006.
 * Copyright Paul A. Bristow 2007.
 *
 * Use, modification and distribution are subject to the
 * Boost Software License, Version 1.0. (See accompanying file
 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
 * ```
 */
 
 'use strict';
 
 // MODULES //
 
 var evalpoly = require( './../../../../base/tools/evalpoly' );
 var erfcinv = require( './../../../../base/special/erfcinv' );
 var abs = require( './../../../../base/special/abs' );
 var exp = require( './../../../../base/special/exp' );
 var ln = require( './../../../../base/special/ln' );
 var sqrt = require( './../../../../base/special/sqrt' );
 var sin = require( './../../../../base/special/sin' );
 var cos = require( './../../../../base/special/cos' );
 var temmeRootFinder = require( './root_finder.js');
 var newtonRaphsonIterate = require( './newton_raphson.js' );
 var polyval1 = require( './polyval_co1.js' );
 var polyval2 = require( './polyval_co2.js' );
 var polyval3 = require( './polyval_co3.js' );
 var polyval4 = require( './polyval_co4.js' );
 var polyval5 = require( './polyval_co5.js' );
 var polyval6 = require( './polyval_co6.js' );
 var polyval7 = require( './polyval_co7.js' );
 var polyval8 = require( './polyval_co8.js' );
 var polyval9 = require( './polyval_co9.js' );
 var polyval10 = require( './polyval_co10.js' );
 var polyval11 = require( './polyval_co11.js' );
 var polyval12 = require( './polyval_co12.js' );
 var polyval13 = require( './polyval_co13.js' );
 
 
 // VARIABLES //
 
 // Workspaces for polynomial coefficients:
 var workspace = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ];
 var terms = [ 0.0, 0.0, 0.0, 0.0 ];
 
 
 // MAIN //
 
 /**
 * Carries out the second method by Temme (described in section 3).
 *
 * ## References
 *
 * -   Temme, N. M. 1992. "Incomplete Laplace Integrals: Uniform Asymptotic Expansion with Application to the Incomplete Beta Function." _Journal of Computational and Applied Mathematics_ 41 (1–2): 1638–63. doi:[10.1016/0377-0427(92)90244-R](https://doi.org/10.1016/0377-0427(92)90244-R).
 *
 * @private
 * @param {number} z - function parameter
 * @param {number} r - function parameter
 * @param {number} theta - function parameter
 * @returns {number} function value
 */
 function temme2( z, r, theta ) {
 	var upper;
 	var lower;
 	var alpha;
 	var roots;
 	var eta0;
 	var eta;
 	var sc7;
 	var sc6;
 	var sc5;
 	var sc4;
 	var sc3;
 	var sc2;
 	var sc;
 	var lu;
 	var s2;
 	var c2;
 	var c;
 	var s;
 	var u;
 	var x;
 
 	// Get first estimate for eta, see Eq 3.9 and 3.10, but note there is a typo in Eq 3.10:
 	eta0 = erfcinv( 2.0*z ) / (-sqrt( r/2.0 ));
 
 	s = sin( theta );
 	c = cos( theta );
 
 	// Now we need to perturb eta0 to get eta, which we do by evaluating the polynomial in 1/r at the bottom of page 151, to do this we first need the error terms e1, e2 e3 which we'll fill into the array "terms".  Since these terms are themselves polynomials, we'll need another array "workspace" to calculate those...
 	terms[ 0 ] = eta0;
 
 	// Some powers of sin(theta) cos(theta) that we'll need later:
 	s2 = s * s;
 	c2 = c * c;
 	sc = s * c;
 	sc2 = sc * sc;
 	sc3 = sc2 * sc;
 	sc4 = sc2 * sc2;
 	sc5 = sc2 * sc3;
 	sc6 = sc3 * sc3;
 	sc7 = sc4 * sc3;
 
 	// Calculate e1 and put it in terms[1], see the middle of page 151:
 	workspace[ 0 ] = ((2.0*s2) - 1.0) / ( 3.0*sc );
 	workspace[ 1 ] = -polyval1( s2 ) / (36.0*sc2);
 	workspace[ 2 ] = polyval2( s2 ) / (1620.0*sc3);
 	workspace[ 3 ] = polyval3( s2 ) / (6480.0*sc4);
 	workspace[ 4 ] = polyval4( s2 ) / (90720.0*sc5);
 	workspace[ 5 ] = 0.0;
 	terms[ 1 ] = evalpoly( workspace, eta0 );
 
 	// Now evaluate e2 and put it in terms[2]:
 	workspace[ 0 ] = -polyval5( s2 ) / (405.0*sc3);
 	workspace[ 1 ] = polyval6( s2 ) / (2592.0*sc4);
 	workspace[ 2 ] = -polyval7( s2 ) / (204120.0*sc5);
 	workspace[ 3 ] = -polyval8( s2 ) / (2099520.0*sc6);
 	workspace[ 4 ] = 0.0;
 	workspace[ 5 ] = 0.0;
 	terms[ 2 ] = evalpoly( workspace, eta0 );
 
 	// And e3, and put it in terms[3]:
 	workspace[ 0 ] = polyval9( s2 ) / (102060.0*sc5);
 	workspace[ 1 ] = -polyval10( s2 ) / (20995200.0*sc6);
 	workspace[ 2 ] = polyval11( s2 ) / (36741600.0*sc7);
 	workspace[ 3 ] = 0.0;
 	workspace[ 4 ] = 0.0;
 	workspace[ 5 ] = 0.0;
 	terms[ 3 ] = evalpoly( workspace, eta0 );
 
 	// Bring the correction terms together to evaluate eta; this is the last equation on page 151:
 	eta = evalpoly( terms, 1.0/r );
 
 	// Now that we have eta we need to back solve for x, we seek the value of x that gives eta in Eq 3.2. The two methods used are described in section 5. Begin by defining a few variables we'll need later:
 	alpha = c / s;
 	alpha *= alpha;
 	lu = ( -( eta*eta )/( 2.0*s2 ) ) + ln(s2) + ( c2*ln(c2)/s2 );
 
 	// Temme doesn't specify what value to switch on here, but this seems to work pretty well:
 	if ( abs(eta) < 0.7 ) {
 		// Small eta use the expansion Temme gives in the second equation of section 5, it's a polynomial in eta:
 		workspace[ 0 ] = s2;
 		workspace[ 1 ] = sc;
 		workspace[ 2 ] = (1.0-(2.0*s2)) / 3.0;
 		workspace[ 3 ] = polyval12( s2 ) / ( 36.0*sc );
 		workspace[ 4 ] = polyval13( s2 ) / ( 270.0*sc2 );
 		workspace[ 5 ] = 0.0;
 		x = evalpoly( workspace, eta );
 	} else {
 		// If eta is large we need to solve Eq 3.2 more directly, begin by getting an initial approximation for x from the last equation on page 155, this is a polynomial in u:
 		u = exp( lu );
 		workspace[ 0 ] = u;
 		workspace[ 1 ] = alpha;
 		workspace[ 2 ] = 0.0;
 		workspace[ 3 ] = 3.0 * alpha * ((3.0*alpha)+1.0) / 6.0;
 		workspace[ 4 ] = 4.0 * alpha * ((4.0*alpha)+1.0) * ((4.0*alpha)+2.0) / 24.0; // eslint-disable-line max-len
 		workspace[ 5 ] = 5.0 * alpha * ((5.0*alpha)+1.0) * ((5.0*alpha)+2.0) * ((5.0*alpha)+3.0) / 120.0; // eslint-disable-line max-len
 		x = evalpoly( workspace, u );
 
 		// At this point we may or may not have the right answer, Eq-3.2 has two solutions for x for any given eta, however the mapping in 3.2 is 1:1 with the sign of eta and x-sin^2(theta) being the same. So we can check if we have the right root of 3.2, and if not switch x for 1-x.  This transformation is motivated by the fact that the distribution is *almost* symmetric so 1-x will be in the right ball park for the solution:
 		if ( (x-s2)*eta < 0.0 ) {
 			x = 1.0 - x;
 		}
 	}
 	// The final step is a few Newton-Raphson iterations to clean up our approximation for x, this is pretty cheap in general, and very cheap compared to an incomplete beta evaluation. The limits set on x come from the observation that the sign of eta and x-sin^2(theta) are the same.
 	if ( eta < 0.0 ) {
 		lower = 0.0;
 		upper = s2;
 	} else {
 		lower = s2;
 		upper = 1.0;
 	}
 	// If our initial approximation is out of bounds then bisect:
 	if ( x < lower || x > upper ) {
 		x = (lower+upper) / 2.0;
 	}
 	roots = temmeRootFinder( -lu, alpha );
 
 	// And iterate:
 	x = newtonRaphsonIterate( roots, x, lower, upper, 32, 100 );
 	return x;
 }
 
 
 // EXPORTS //
 
 module.exports = temme2;