import * as bfjs from "./bf.js";
/**
* High-precision Root Finding using the Brent-Dekker Method (similar to MATLAB's fzero).
*
* This algorithm combines the reliability of Bisection, the speed of the Secant method,
* and the high-order convergence of Inverse Quadratic Interpolation (IQI).
* It guarantees global convergence while achieving superlinear convergence rates near the root.
*
* @param {Function} f - The target function to find the root of.
* Must accept a BigFloat argument and return a BigFloat result.
* The function values at the endpoints must have opposite signs (f(_a) * f(_b) <= 0).
*
* @param {number|string|BigFloat} _a - The start of the search interval (or first initial guess).
* @param {number|string|BigFloat} _b - The end of the search interval (or second initial guess).
*
*
* @param {Object} [info={}] - Configuration and Status object.
* This object configures the execution parameters and is updated in-place
* with statistical data during and after execution.
* @param {number|string|BigFloat} [info._e=1e-30] - Absolute Error Tolerance.
* Convergence is considered achieved when the interval width or step size falls below this value.
*
* @param {number|string|BigFloat} [info._re=info._e] - Relative Error Tolerance.
* Used to handle convergence for large values. Defaults to the absolute tolerance.
* The effective tolerance is calculated as: `tol = |b| * _re + _e`.
*
* // --- Input Configuration Properties ---
* @param {number} [info.max_step=200] - Maximum number of iterations allowed.
* If this limit is reached without convergence, the function returns null and logs a warning.
* @param {number} [info.max_time=60000] - Maximum execution time in milliseconds.
* Prevents the function from hanging in infinite loops or extremely slow computations.
* @param {Function} [info.cb] - Optional callback function.
* If defined, this function is called after every iteration. Useful for updating UI progress or logging.
* @param {boolean} [info.debug] - Optional flag to enable debug logging (implementation dependent).
*
* // --- Output Status Properties (Updated during execution) ---
* @param {BigFloat|null} info.result - The final found root.
* Returns a BigFloat if converged, or null if failed.
* @param {BigFloat} info.lastresult - The result of the last iteration (current best guess `b`).
* Even if convergence fails, this contains the value closest to the root when execution stopped.
* @param {string} info.eff_result - String representation of the result based on effective precision.
* Generated by truncating `lastresult` according to `eff_decimal_precision`.
* @param {number} info.steps - The number of iterations currently executed.
* @param {number} info.exectime - The elapsed execution time in milliseconds.
* @param {BigFloat} info.error - The estimated error bound.
* Typically represents half the width of the current search interval (`xm`).
* @param {BigFloat} info.residual - The absolute value of the function at the current best guess: `|f(x)|`.
* Ideally, this value should be close to zero.
* @param {number} info.eff_decimal_precision - Estimated number of significant decimal digits.
* Calculated as `-log10(error)`.
* @param {Function} info.toString - Helper method.
* Returns a formatted string containing steps, error, residual, and execution time.
*
* @returns {BigFloat|null}
* Returns the BigFloat root if the tolerance criteria are met.
* Returns `null` if the maximum steps or time limit is exceeded (a warning is logged to the console).
*/
export function fzero(f, _a, _b, info = {}) {
let _e = info._e ?? 1e-30;
let _re = info._re ?? _e;
// 1. Configuration & Constants
// Default steps for root finding (usually faster than integration, but high precision needs more)
let max_step = info.max_step || 200,
max_time = info.max_time || 60000;
if (typeof(_e) !== 'number' || typeof(_re) !== 'number' || typeof(info) !== "object") {
throw new Error("arguments error");
}
const start_time = new Date().getTime();
// Cache BigFloat constants to avoid repeated object creation
const bf_zero = bfjs.zero;
const bf_one = bfjs.one;
const bf_two = bfjs.two;
const bf_three = bfjs.three;
const bf_half = bfjs.half;
// Initial Interval
let a = bfjs.bf(_a);
let b = bfjs.bf(_b);
let fa = f(a);
let fb = f(b);
// Initial check: f(a) and f(b) must have opposite signs
if (fa.mul(fb).cmp(bf_zero) > 0) {
throw new Error("Function values at the interval endpoints must differ in sign.");
}
// `c` is the "contrapoint" (last point with opposite sign to b).
// In Brent's method, the root is always bracketed between b and c.
let c = a;
let fc = fa;
let d = b.sub(a); // Step size
let e = d; // Previous step size
// Tolerances
let tol_act = bfjs.bf(_e);
let tol_rel = bfjs.bf(_re);
// Helper to sync local variables to info object (Performance optimization)
const updateInfo = (iter, errorBound, residual) => {
info.exectime = new Date().getTime() - start_time;
info.steps = iter;
info.lastresult = b; // Current best guess
info.residual = residual;
info.error = errorBound;
// Calculate effective decimal precision
// Handle log(0) explicitly
if (errorBound.cmp(bf_zero) === 0) {
info.eff_decimal_precision = bfjs.decimalPrecision();
} else {
info.eff_decimal_precision = Math.floor(-errorBound.log().f64() / Math.log(10));
}
if (info.eff_decimal_precision <= 0) {
info.eff_decimal_precision = 0;
info.eff_result = '';
} else {
let limit = bfjs.decimalPrecision();
let prec = info.eff_decimal_precision > limit ? limit : info.eff_decimal_precision;
info.eff_result = b.toString(10, prec);
}
};
// Format output string
info.toString = function() {
return `root=${this.eff_result},
residual=${this.residual ? this.residual.toString(10, 3) : 'N/A'},
steps=${this.steps}/${max_step},
error=${this.error ? this.error.toString(10, 3) : 'N/A'},
eff_prec=${this.eff_decimal_precision},
exectime=${this.exectime}/${max_time}`;
};
// 2. Main Loop
for (let iter = 1; iter <= max_step; ++iter) {
// Ensure |f(b)| <= |f(c)|. b is always the best guess.
if (fb.abs().cmp(fc.abs()) > 0) {
a = c; fa = fc;
c = b; fc = fb;
b = a; fb = fa;
}
// Convergence tolerances
// tol1 = |b| * tol_rel + tol_act
let tol1 = b.abs().mul(tol_rel).add(tol_act);
let xm = c.sub(b).mul(bf_half); // Midpoint relative to b
let xm_abs = xm.abs();
// Check Convergence
// 1. Interval size is smaller than tolerance
// 2. Function value is exactly zero
if (xm_abs.cmp(tol1) <= 0 || fb.cmp(bf_zero) === 0) {
updateInfo(iter, xm_abs, fb.abs());
info.result = b;
return b;
}
// Call callback if exists (sync info first)
if (info.cb) {
updateInfo(iter, xm_abs, fb.abs());
info.cb();
}
// Check timeout
if (new Date().getTime() - start_time > max_time) {
updateInfo(iter, xm_abs, fb.abs());
console.log("fzero: Timeout reached.");
info.result = null;
return null;
}
// Determine Step Strategy
// Attempt Inverse Quadratic Interpolation (IQI) or Secant if:
// 1. Previous step was large enough (|e| >= tol1)
// 2. Function value decreased enough (|fa| > |fb|)
if (e.abs().cmp(tol1) >= 0 && fa.abs().cmp(fb.abs()) > 0) {
let s = fb.div(fa);
let p, q;
// If a == c, we only have 2 distinct points -> Secant Method
if (a.cmp(c) === 0) {
// p = 2 * xm * s
// q = 1 - s
p = xm.mul(bf_two).mul(s);
q = bf_one.sub(s);
} else {
// We have 3 distinct points -> Inverse Quadratic Interpolation
let q_iqi = fa.div(fc);
let r_iqi = fb.div(fc);
// p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0))
let term1 = xm.mul(bf_two).mul(q_iqi).mul(q_iqi.sub(r_iqi));
let term2 = b.sub(a).mul(r_iqi.sub(bf_one));
p = s.mul(term1.sub(term2));
// q = (q - 1.0) * (r - 1.0) * (s - 1.0)
q = q_iqi.sub(bf_one).mul(r_iqi.sub(bf_one)).mul(s.sub(bf_one));
}
// Adjust signs so q > 0 to simplify division logic
if (p.cmp(bf_zero) > 0) {
q = q.neg();
} else {
p = p.neg();
}
q = q.abs();
// Validate Interpolation
// Condition 1: Must be within the interval (bounded by 3*xm*q - |tol1*q|)
// Condition 2: Must converge faster than bisection (p < |0.5*e*q|)
let cond1_bound = xm.mul(bf_three).mul(q).sub(tol1.mul(q).abs());
let cond2_bound = e.mul(q).abs().mul(bf_half);
if (p.mul(bf_two).cmp(cond1_bound) < 0 && p.cmp(cond2_bound) < 0) {
// Accept Interpolation Step
e = d;
d = p.div(q);
} else {
// Reject, use Bisection
d = xm;
e = d;
}
} else {
// Bounds too slow or not converging fast enough -> Bisection
d = xm;
e = d;
}
// Apply Step
a = b;
fa = fb;
// Ensure step is at least tol1 (to avoid underflow/stagnation near root)
if (d.abs().cmp(tol1) > 0) {
b = b.add(d);
} else {
// b += sign(xm) * tol1
let signXm = xm.cmp(bf_zero) >= 0 ? bf_one : bf_one.neg();
b = b.add(tol1.mul(signXm));
}
fb = f(b);
// Maintain Bracket: If fb and fc have same sign, reset c to a
// (Root must always be between b and c)
if (fb.cmp(bf_zero) === 0 || fb.mul(fc).cmp(bf_zero) > 0) {
c = a;
fc = fa;
d = b.sub(a);
e = d;
}
}
// 3. Exit (Max steps reached)
// Update info one last time before failure
let final_xm = c.sub(b).mul(bf_half).abs();
updateInfo(max_step, final_xm, fb.abs());
console.log(`fzero: Failed to converge after ${max_step} steps. Residual: ${fb.toString(10,3)}`);
info.result = null;
return null;
};