/**
 * Searches the interval from <tt>lowerLimit</tt> to <tt>upperLimit</tt>
 * for a root (i.e., zero) of the function <tt>func</tt> with respect to
 * its first argument using Brent's method root-finding algorithm.
 *
 * Translated from zeroin.c in http://www.netlib.org/c/brent.shar.
 *
 * Copyright (c) 2012 Borgar Thorsteinsson <borgar@borgar.net>
 * MIT License, http://www.opensource.org/licenses/mit-license.php
 *
 * @param {function} function for which the root is sought.
 * @param {number} the lower point of the interval to be searched.
 * @param {number} the upper point of the interval to be searched.
 * @param {number} the desired accuracy (convergence tolerance).
 * @param {number} the maximum number of iterations.
 * @returns an estimate for the root within accuracy.
 *
 */
export default function uniroot(
	func,
	lowerLimit,
	upperLimit,
	errorTol,
	maxIter
) {
	var a = lowerLimit,
		b = upperLimit,
		c = a,
		fa = func(a),
		fb = func(b),
		fc = fa,
		tol_act, // Actual tolerance
		new_step, // Step at this iteration
		prev_step, // Distance from the last but one to the last approximation
		p, // Interpolation step is calculated in the form p/q; division is delayed until the last moment
		q

	errorTol = errorTol || 0
	maxIter = maxIter || 1000

	while (maxIter-- > 0) {
		prev_step = b - a

		if (Math.abs(fc) < Math.abs(fb)) {
			// Swap data for b to be the best approximation
			;(a = b), (b = c), (c = a)
			;(fa = fb), (fb = fc), (fc = fa)
		}

		tol_act = 1e-15 * Math.abs(b) + errorTol / 2
		new_step = (c - b) / 2

		if (Math.abs(new_step) <= tol_act || fb === 0) {
			return b // Acceptable approx. is found
		}

		// Decide if the interpolation can be tried
		if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
			// If prev_step was large enough and was in true direction, Interpolatiom may be tried
			var t1, cb, t2
			cb = c - b
			if (a === c) {
				// If we have only two distinct points linear interpolation can only be applied
				t1 = fb / fa
				p = cb * t1
				q = 1.0 - t1
			} else {
				// Quadric inverse interpolation
				;(q = fa / fc), (t1 = fb / fc), (t2 = fb / fa)
				p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1))
				q = (q - 1) * (t1 - 1) * (t2 - 1)
			}

			if (p > 0) {
				q = -q // p was calculated with the opposite sign; make p positive
			} else {
				p = -p // and assign possible minus to q
			}

			if (
				p < 0.75 * cb * q - Math.abs(tol_act * q) / 2 &&
				p < Math.abs(prev_step * q / 2)
			) {
				// If (b + p / q) falls in [b,c] and isn't too large it is accepted
				new_step = p / q
			}

			// If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
		}

		if (Math.abs(new_step) < tol_act) {
			// Adjust the step to be not less than tolerance
			new_step = new_step > 0 ? tol_act : -tol_act
		}

		;(a = b), (fa = fb) // Save the previous approx.
		;(b += new_step), (fb = func(b)) // Do step to a new approxim.

		if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
			;(c = a), (fc = fa) // Adjust c for it to have a sign opposite to that of b
		}
	}
}

/*
var test_counter;
function f1 (x) { test_counter++; return (Math.pow(x,2)-1)*x - 5; }
function f2 (x) { test_counter++; return Math.cos(x)-x; }
function f3 (x) { test_counter++; return Math.sin(x)-x; }
function f4 (x) { test_counter++; return (x + 3) * Math.pow(x - 1, 2); }
[
  [f1, 2, 3],
  [f2, 2, 3],
  [f2, -1, 3],
  [f3, -1, 3],
  [f4, -4, 4/3]
].forEach(function (args) {
  test_counter = 0;
  var root = uniroot.apply( pv, args );
  ;;;console.log( 'uniroot:', args.slice(1), root, test_counter );
})
*/