Function optimal_neumann_terms 

pub fn optimal_neumann_terms(
    coherence: Precision,
    b_inf_norm: Precision,
    min_diag: Precision,
    tolerance: Precision,
) -> Option<usize>

Minimum number of Neumann terms required to hit tolerance on a single-entry solve, given the matrix’s coherence margin and the magnitudes of the RHS and any perturbation.

Math

For strictly DD A = D - O with coherence c, the Neumann series A⁻¹ b = Σ_k y_k satisfies ‖y_k‖_∞ ≤ ρ^k · ‖y_0‖_∞ where the spectral radius of D⁻¹O is bounded by ρ ≤ 1 - c. To get per-entry error below tolerance we need

  k ≥ log(‖y_0‖_∞ / tolerance) / log(1 / ρ)
      = log(‖b‖_∞ / (min_diag · tolerance)) / log(1 / (1 - c))

This function picks the smallest such k (rounded up). Saturates at 64 so a callers-provided pathological pair can’t blow up the iteration count.

Returns

• Some(k) — recommended Neumann term count, clamped to [1, 64].
• None — non-strict-DD input (coherence <= 0), or non-positive tolerance, or non-positive min_diag. Caller should fall back to a hand-picked max_terms or to a full solve.

Edge cases

• b_inf_norm == 0 → returns Some(1) (we still want at least one term to confirm the zero answer).
• Inputs that produce k > 64 are clamped (the bound is conservative; real matrices rarely need that many).

Examples

// For coherence 0.5, ‖b‖ = 10, min_diag = 5, tolerance = 1e-8:
//   k ≥ log(10 / (5 · 1e-8)) / log(2) = log2(2e8) ≈ 27.6 → k = 28
let k = optimal_neumann_terms(
    /*coherence=*/ 0.5,
    /*b_inf_norm=*/ 10.0,
    /*min_diag=*/   5.0,
    /*tolerance=*/  1e-8,
).unwrap();
assert!(k >= 28 && k <= 30);