fix: avoid u64 overflow in solve_itp for large nmax values

[greysquirr3l]

Problem

When epsilon is very small (e.g. 1e-24), nmax in solve_itp_fallible can exceed 63. The current expression:

let mut scaled_epsilon = epsilon * (1u64 << nmax) as f64;

overflows u64 when nmax >= 64:

• Debug builds: panic with "attempt to shift left with overflow"
• Release builds: silent wrong result — x86 masks shift amount to 6 bits, so << 64 wraps to << 0, producing scaled_epsilon = epsilon * 1 instead of epsilon * 2^nmax

How nmax gets large

let n1_2 = (((b - a) / epsilon).log2().ceil() - 1.0).max(0.0) as usize;
let nmax = n0 + n1_2;

For b - a = 1.0 and epsilon = 1e-24, log2(1e24) ≈ 79.7, giving nmax ≈ 79.

Fix

Replace the integer shift with f64::exp2():

let mut scaled_epsilon = epsilon * (nmax as f64).exp2();

This is mathematically identical (2^nmax) but uses floating-point arithmetic. f64 handles nmax up to ~1023 before saturating to +inf — well beyond any epsilon that makes physical sense.

Test

Added test_solve_itp_large_nmax_does_not_overflow with epsilon = 1e-24, which forces nmax ≈ 79. This panics on the original code in debug mode and returns the wrong root in release mode.

[greysquirr3l]

Ope, investigating CI issue.

[Copilot]

Pull request overview

Fixes an overflow bug in the ITP root solver when epsilon is extremely small and the derived nmax exceeds 63, preventing debug-mode panics and release-mode misbehavior.

Changes:

• Replaces a u64 left-shift (1u64 << nmax) with a floating-point exp2 computation when scaling epsilon.
• Adds a regression test that exercises a case where nmax becomes large (e.g. epsilon = 1e-24).

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[mlwilkerson]

Oh, I think this may be the fix for issue #602 I submitted recently.

[raphlinus]

I'm going to dig deeply into this set of problems soon. My immediate question is whether nmax > 63 is ever valid. A back-of-the-envelope analysis suggests that with bisection you're never going to get more than the number of mantissa bits (52 or 53 bits), so I'm thinking of just clamping nmax to that. That may well mean not being able to satisfy the precision request implied by epsilon.

Clearly the existing code is wrong and needs a fix. I just need to find out where these too-small epsilons are coming from and what the best behavior for the solver should be in those cases.