time-to-botec

halley_iterate.js (6068B)

---

 /**
 * @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_62_0/boost/math/tools/roots.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * Copyright John Maddock 2006.
 *
 * 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 abs = require( './../../../../base/special/abs' );
 var ldexp = require( './../../../../base/special/ldexp' );
 var sign = require( './../../../../base/special/signum' );
 var max = require( './../../../../base/special/max' );
 var MAX_VALUE = require( '@stdlib/constants/float64/max' );
 
 
 // MAIN //
 
 /**
 * Performs root finding via third order Halley iteration.
 *
 * @private
 * @param {Array} fun - array of function and its first two derivatives
 * @param {number} guess - initial starting value
 * @param {number} minimum - minimum possible value for the result, used as initial lower bracket
 * @param {number} maximum - maximum possible value for the result, used as initial upper bracket
 * @param {PositiveInteger} digits - desired number of binary digits
 * @param {PositiveInteger} maxIter - maximum number of iterations
 * @returns {number} function value
 */
 function halleyIterate( fun, guess, minimum, maximum, digits, maxIter ) {
 	var convergence;
 	var outOfBounds;
 	var delta1;
 	var delta2;
 	var factor;
 	var result;
 	var f0Last;
 	var count;
 	var delta;
 	var denom;
 	var diff;
 	var num;
 	var res;
 	var f0;
 	var f1;
 	var f2;
 
 	f0 = 0.0;
 	outOfBounds = false;
 	result = guess;
 	factor = ldexp( 1.0, 1.0-digits );
 	delta = max( 10000000*guess, 10000000 );  // Arbitrarily large delta...
 	f0Last = 0;
 	delta1 = delta;
 	delta2 = delta;
 
 	count = maxIter;
 	do {
 		f0Last = f0;
 		delta2 = delta1;
 		delta1 = delta;
 		res = fun( result);
 		f0 = res[ 0 ];
 		f1 = res[ 1 ];
 		f2 = res[ 2 ];
 		count -= 1;
 
 		if ( f0 === 0.0 ) {
 			break;
 		}
 		if ( f1 === 0.0 ) {
 			// Oops zero derivative!!!
 			if ( f0Last === 0.0 ) {
 				// Must be the first iteration, pretend that we had a previous one at either min or max:
 				if ( result === minimum ) {
 					guess = maximum;
 				} else {
 					guess = minimum;
 				}
 				f0Last = fun( guess );
 				delta = guess - result;
 			}
 			if ( sign( f0Last ) * sign( f0 ) < 0 ) {
 				// We've crossed over so move in opposite direction to last step:
 				if ( delta < 0 ) {
 					delta = ( result-minimum ) / 2.0;
 				} else {
 					delta = ( result-maximum ) / 2.0;
 				}
 			// Move in same direction as last step:
 			} else if ( delta < 0 ) {
 				delta = (result-maximum) / 2.0;
 			} else {
 				delta = (result-minimum) / 2.0;
 			}
 		} else if ( f2 === 0.0 ) {
 			delta = f0 / f1;
 		} else {
 			denom = 2.0 * f0;
 			num = ( 2.0 * f1 ) - ( f0 * ( f2 / f1 ) );
 			if ( abs(num) < 1.0 && ( abs(denom) >= abs(num) * MAX_VALUE ) ) {
 				// Possible overflow, use Newton step:
 				delta = f0 / f1;
 			} else {
 				delta = denom / num;
 			}
 			if ( delta * f1 / f0 < 0.0 ) {
 				// Probably cancellation error, try a Newton step instead:
 				delta = f0 / f1;
 				if ( abs(delta) > 2.0 * abs(guess) ) {
 					delta = ( (delta < 0.0) ? -1.0 : 1.0 ) * 2.0 * abs( guess );
 				}
 			}
 		}
 		convergence = abs( delta / delta2 );
 		if ( convergence > 0.8 && convergence < 2.0 ) {
 			// Last two steps haven't converged, try bisection:
 			delta = ( delta > 0.0 ) ? ( result-minimum )/2.0 : ( result-maximum )/2.0; // eslint-disable-line max-len
 			if ( abs(delta) > result ) {
 				delta = sign( delta ) * result; // Protect against huge jumps!
 			}
 			// Reset delta2 so that this branch will *not* be taken on the next iteration:
 			delta2 = delta * 3.0;
 		}
 		guess = result;
 		result -= delta;
 
 		// Check for out of bounds step:
 		if ( result < minimum ) {
 			if (
 				abs(minimum) < 1 &&
 				abs(result) > 1 &&
 				( MAX_VALUE / abs(result) < abs(minimum) )
 			) {
 				diff = 1000.0;
 			} else {
 				diff = result / minimum;
 			}
 			if ( abs(diff) < 1.0 ) {
 				diff = 1.0 / diff;
 			}
 			if ( !outOfBounds && diff > 0.0 && diff < 3.0 ) {
 				// Only a small out of bounds step, let's assume that the result is probably approximately at minimum:
 				delta = 0.99 * (guess - minimum);
 				result = guess - delta;
 				outOfBounds = true; // Only take this branch once!
 			} else {
 				delta = (guess - minimum) / 2.0;
 				result = guess - delta;
 				if ( result === minimum || result === maximum ) {
 					break;
 				}
 			}
 		} else if ( result > maximum ) {
 			if (
 				abs(maximum) < 1.0 &&
 				abs(result) > 1.0 &&
 				MAX_VALUE / abs(result) < abs(maximum)
 			) {
 				diff = 1000.0;
 			} else {
 				diff = result / maximum;
 			}
 			if ( abs(diff) < 1.0 ) {
 				diff = 1.0 / diff;
 			}
 			if ( !outOfBounds && diff > 0.0 && diff < 3.0 ) {
 				// Only a small out of bounds step, let's assume that the result is probably approximately at minimum:
 				delta = 0.99 * (guess - maximum);
 				result = guess - delta;
 				outOfBounds = true; // Only take this branch once!
 			} else {
 				delta = ( guess - maximum ) / 2.0;
 				result = guess - delta;
 				if ( result === minimum || result === maximum ) {
 					break;
 				}
 			}
 		}
 		// Update brackets:
 		if ( delta > 0.0 ) {
 			maximum = guess;
 		} else {
 			minimum = guess;
 		}
 	} while ( count && ( abs(result * factor) < abs(delta) ) );
 
 	return result;
 }
 
 
 // EXPORTS //
 
 module.exports = halleyIterate;