/**
* @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;