export class NoBracketError extends Error {}

/**
 * Widens [a, b] outward geometrically until f changes sign.
 * The forward is monotonic in every solvable input, so one sign change ⇒ one root.
 */
export function expandBracket(
  f: (x: number) => number,
  a: number,
  b: number,
  maxIter = 50,
  factor = 1.6,
): [number, number] {
  let fa = f(a);
  let fb = f(b);
  for (let i = 0; i < maxIter; i++) {
    if (Number.isFinite(fa) && Number.isFinite(fb) && fa * fb <= 0) return [a, b];
    if (Math.abs(fa) < Math.abs(fb)) {
      a += factor * (a - b);
      fa = f(a);
    } else {
      b += factor * (b - a);
      fb = f(b);
    }
  }
  throw new NoBracketError(`No sign change found in [${a}, ${b}]`);
}

/** Brent's method (inverse quadratic interpolation + secant + bisection fallback). */
export function brent(
  f: (x: number) => number,
  lo: number,
  hi: number,
  tol = 1e-12,
  maxIter = 200,
): number {
  let a = lo;
  let b = hi;
  let fa = f(a);
  let fb = f(b);
  if (fa === 0) return a;
  if (fb === 0) return b;
  if (fa * fb > 0) throw new NoBracketError(`Root not bracketed in [${lo}, ${hi}]`);

  let c = a;
  let fc = fa;
  let d = b - a;
  let e = d;

  for (let i = 0; i < maxIter; i++) {
    if (fb * fc > 0) {
      c = a;
      fc = fa;
      d = e = b - a;
    }
    if (Math.abs(fc) < Math.abs(fb)) {
      a = b; b = c; c = a;
      fa = fb; fb = fc; fc = fa;
    }
    const tol1 = 2 * Number.EPSILON * Math.abs(b) + 0.5 * tol;
    const xm = 0.5 * (c - b);
    if (Math.abs(xm) <= tol1 || fb === 0) return b;

    if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
      const s = fb / fa;
      let p: number;
      let q: number;
      if (a === c) {
        p = 2 * xm * s;
        q = 1 - s;
      } else {
        const qq = fa / fc;
        const r = fb / fc;
        p = s * (2 * xm * qq * (qq - r) - (b - a) * (r - 1));
        q = (qq - 1) * (r - 1) * (s - 1);
      }
      if (p > 0) q = -q;
      p = Math.abs(p);
      if (2 * p < Math.min(3 * xm * q - Math.abs(tol1 * q), Math.abs(e * q))) {
        e = d;
        d = p / q;
      } else {
        d = xm;
        e = d;
      }
    } else {
      d = xm;
      e = d;
    }
    a = b;
    fa = fb;
    b += Math.abs(d) > tol1 ? d : Math.sign(xm) * tol1;
    fb = f(b);
  }
  return b;
}