Function fit_pow3 

pub fn fit_pow3(points: &[(f64, f64)]) -> Option<Pow3Fit>

Fits the Domhan pow-3 law y(t) = c − a·t^(−α) to points (each (t, y) with t ≥ 1) by in-house Levenberg–Marquardt least squares.

The residual of point i is r_i = y_i − (c − a·t_i^(−α)), and the fit minimises Σ r_i². Writing p_i = t_i^(−α), the residual Jacobian columns are

∂r/∂a = p_i,   ∂r/∂α = −a·p_i·ln(t_i),   ∂r/∂c = −1

so each iteration forms the 3×3 Gauss–Newton normal matrix H = JᵀJ and gradient g = Jᵀr, solves the damped system (H + λ·diag H) δ = −g by Cramer’s rule (solve3), and accepts the step when it lowers the sum of squares (shrinking λ) or rejects it (growing λ) — plain Marquardt damping. It is initialised with α₀ = 1, c₀ just above the largest observed value (the asymptote sits above a rising curve), and a₀ chosen to match the first point exactly.

Returns None for a degenerate, failed, or unsupported fit — fewer than three points, a non-finite input, a singular normal matrix, a non-finite result, or an unsupported extrapolation (a near-linear curve with decay exponent below MIN_ALPHA, or an asymptote more than REMAINING_FACTOR× the observed range above the data) — so a caller that cannot model a curve, or that can only mirage one, declines to act on it rather than inventing a number This is the fail-safe that keeps a still rising trial’s inflated projection from freezing a genuinely-better plateaued one. a and α are kept strictly positive throughout.