time-to-botec

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

newton_raphson.js (3856B)

      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_62_0/boost/math/tools/roots.hpp}. The implementation has been modified for JavaScript.
     22 *
     23 * ```text
     24 * 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 sign = require( './../../../../base/special/signum' );
     37 var abs = require( './../../../../base/special/abs' );
     38 var ldexp = require( './../../../../base/special/ldexp' );
     39 var MAX_VALUE = require( '@stdlib/constants/float64/max' );
     40 
     41 
     42 // MAIN //
     43 
     44 /**
     45 * Performs root finding via second order Newton-Raphson iteration.
     46 *
     47 * @private
     48 * @param {Array} fun - two-element array of the function and its first derivative
     49 * @param {number} guess - initial starting value.
     50 * @param {number} min - minimum possible value for the result,used as initial lower bracket.
     51 * @param {number} max - maximum possible value for the result, used as initial upper bracket.
     52 * @param {PositiveInteger} digits - desired number of binary digits
     53 * @param {PositiveInteger} maxIter - maximum number of iterations
     54 * @returns {number} function value
     55 */
     56 function newtonRaphsonIterate( fun, guess, min, max, digits, maxIter ) {
     57 	var f0last;
     58 	var delta1;
     59 	var delta2;
     60 	var factor;
     61 	var result;
     62 	var count;
     63 	var delta;
     64 	var res;
     65 	var f0;
     66 	var f1;
     67 
     68 	f0 = 0.0;
     69 	f0last = 0.0;
     70 	result = guess;
     71 
     72 	factor = ldexp( 1.0, 1.0 - digits );
     73 	delta = MAX_VALUE;
     74 	delta1 = MAX_VALUE;
     75 	delta2 = MAX_VALUE;
     76 
     77 	count = maxIter;
     78 	do {
     79 		f0last = f0;
     80 		delta2 = delta1;
     81 		delta1 = delta;
     82 		res = fun(result);
     83 		f0 = res[ 0 ];
     84 		f1 = res[ 1 ];
     85 		count -= 1;
     86 		if ( f0 === 0.0 ) {
     87 			break;
     88 		}
     89 		if ( f1 === 0.0 ) {
     90 			// Oops zero derivative!!!
     91 			if ( f0last === 0.0 ) {
     92 				// Must be the first iteration, pretend that we had a previous one at either min or max:
     93 				if ( result === min ) {
     94 					guess = max;
     95 				} else {
     96 					guess = min;
     97 				}
     98 				f0last = fun( guess );
     99 				delta = guess - result;
    100 			}
    101 			if ( sign(f0last) * sign(f0) < 0 ) {
    102 				// We've crossed over so move in opposite direction to last step:
    103 				if ( delta < 0 ) {
    104 					delta = (result - min) / 2.0;
    105 				} else {
    106 					delta = (result - max) / 2.0;
    107 				}
    108 			} else if ( delta < 0 ) {
    109 				delta = (result - max) / 2.0;
    110 			} else {
    111 				delta = (result - min) / 2.0;
    112 			}
    113 		} else {
    114 			delta = f0 / f1;
    115 		}
    116 		if ( abs(delta * 2.0) > abs(delta2) ) {
    117 			// Last two steps haven't converged, try bisection:
    118 			delta = ( delta > 0.0 ) ? (result-min) / 2.0 : (result-max) / 2.0;
    119 		}
    120 		guess = result;
    121 		result -= delta;
    122 		if ( result <= min ) {
    123 			delta = 0.5 * (guess - min);
    124 			result = guess - delta;
    125 			if ( result === min || result === max ) {
    126 				break;
    127 			}
    128 		} else if ( result >= max ) {
    129 			delta = 0.5 * (guess - max);
    130 			result = guess - delta;
    131 			if ( result === min || result === max ) {
    132 				break;
    133 			}
    134 		}
    135 		// Update brackets:
    136 		if ( delta > 0.0 ) {
    137 			max = guess;
    138 		} else {
    139 			min = guess;
    140 		}
    141 	}
    142 	while ( count && ( abs(result * factor) < abs(delta) ) );
    143 
    144 	return result;
    145 }
    146 
    147 
    148 // EXPORTS //
    149 
    150 module.exports = newtonRaphsonIterate;