Make solve_itp terminate, and build its schedule without an overflow

[kooshi]

AI disclosure: these bugs were diagnosed and resolved by Claude Opus 5, working under my guidance to be exceptionally rigorous. I ran into these issues in a real codebase, needed a fix, and have verified that this resolves the issue.
Everything beyond this point is written by Claude.

---

inv_arclen hangs on long curves and panics on longer ones. Two lines, kurbo 0.13.1, debug build,
and this never comes back:

let c = CubicBez::new((0., 0.), (5e7, 0.), (5e7, 1e7), (1e8, 1e7));
c.inv_arclen(c.arclen(1e-9) * 0.171, 1e-9);

It sits there burning a full core and allocating nothing. The curve is about 1e8 long, and 0.171 is
nothing special: of 999 evenly spaced positions along it, 199 hang the same way. On a curve ten
times longer, 248 of 999. Short curves are fine.

Make the curve longer still and you get a panic instead. In release, where there is no overflow
check, that one hangs too:

let c = CubicBez::new((0., 0.), (5e9, 0.), (5e9, 1e9), (1e10, 1e9));
c.inv_arclen(c.arclen(1e-9) * 0.3137, 1e-9);
// panicked at common.rs:713:40: attempt to shift left with overflow

I ran into this through fit_to_bezpath. Its fitter calls inv_arclen(_, 1e-9) on candidate
curves, one candidate came out long enough, and the fit never returned.

Below is what I think is going on and what I did about it. You know this code better than I do, so
treat the diagnosis as a proposal too.

How the loop works, for anyone who has not read the paper

solve_itp finds where a function crosses zero. You hand it a bracket [a, b] with the function
negative at one end and positive at the other, and it repeatedly picks a point inside, evaluates
there, and throws away whichever half cannot contain the crossing. The bracket closes in. It stops
once the bracket is narrower than 2 * epsilon, so epsilon is how close to the crossing you are
asking it to get.

Bisection would just take the midpoint every time. Safe, and exactly one bit per call. Regula falsi
draws a straight line between the two endpoints and takes where that line crosses zero, which is
usually much faster and occasionally much worse.

ITP is a way of getting the second without giving up the first. Each pass it works out the regula
falsi point xf, nudges it a little toward the midpoint (xt), and then, and this is the part that
matters here, refuses to move it further than r away from the midpoint:

let x1_2 = 0.5 * (a + b);
let r = scaled_epsilon - 0.5 * (b - a);
// ... xf and xt worked out from the endpoints ...
let xitp = if (xt - x1_2).abs() <= r { xt } else { x1_2 - r.copysign(sigma) };

r is the interesting quantity. scaled_epsilon is the widest the bracket is still allowed to be
at this point in the run: it starts at epsilon * 2^nmax, where nmax is roughly how many halvings
it would take to get from the starting bracket down to epsilon, plus a few spare passes (n0),
and it halves every iteration. Subtract the half-width you actually have and you get r, which is
how far ahead of that schedule you are running.

So r is a budget. Ahead of schedule and it is positive, which buys room to move away from the
midpoint and trust the interpolated point. Behind schedule and it goes to zero, xitp is pinned to
the midpoint, and the algorithm turns back into plain bisection. That is the trick: the guessing is
only allowed while there is slack to pay for it.

Bug one: the loop can have no way out

inv_arclen works out the epsilon it needs like this:

let epsilon = accuracy / total_arclen;

Which is the right conversion. You asked for the answer to within accuracy units of arc length,
the answer comes back as a parameter t in [0, 1], so the tolerance on t is accuracy divided
by how long the curve is. On a curve 1e8 long at accuracy = 1e-9, that is about 1e-17.

The trouble is that t is an f64. Around 0.125 the gap between one f64 and the next is about
2.8e-17. So the bracket is being asked to get narrower than the gap between the numbers it is made
of, which it cannot do. It shrinks until a and b are neighbours and then stops, because there is
nothing left in between.

Meanwhile while b - a > 2.0 * epsilon is still true, and stays true forever. And the point the
loop picks now rounds onto a or b, so a = xitp assigns a to itself. Nothing changes, so the
next pass computes the same point, and the one after that. Here is a run of it, from a small program
that prints what the objective gets asked for:

  a = 1.25186695533393788082e-1
  b = 1.25186695533393815838e-1
  b - a       = 2.775558e-17   (one f64 step, nothing can sit between them)
  2 * epsilon = 2.000000e-17   (the loop exits only when b - a drops below this)
...
call  2520000000  x = 1.25186695533393788082e-1  same x 2519999944 times in a row

Two and a half billion calls, all at the same point. It stopped making progress at call 56.

There is a second way in, with no curve involved at all. b - a is never negative and 2.0 * epsilon is always negative for a negative epsilon, so the exit test cannot go false for any
bracket at all. That takes four arguments and a minus sign.

Bug two: 1u64 << nmax runs out of room

nmax never gets bounded at 63:

let n1_2 = (((b - a) / epsilon).log2().ceil() - 1.0).max(0.0) as usize;
let nmax = n0 + n1_2;
let mut scaled_epsilon = epsilon * (1u64 << nmax) as f64;

A curve about 1e10 long gives epsilon = 9.9e-20, and nmax comes out at 64. Debug builds panic.
Release builds have no check, so the shift wraps to nmax % 64 and scaled_epsilon starts out a
factor of 2^64 too small. Then r is negative on every pass, xitp is stuck on the midpoint, and
the call lands in bug one instead. So the release build is the worse of the two, not the milder one.

There is also an overflow one line up. epsilon = 0.0 makes (b - a) / epsilon infinite, inf as usize saturates to usize::MAX, and n0 + n1_2 overflows the addition before the shift is
reached. solve_itp is public and takes epsilon directly, so that one needs no curve either.

The doc comment does ask callers to keep epsilon above 2^-63 times b - a, and bug two is that
line being crossed. Two things about it. The caller crossing it is kurbo, since someone using
inv_arclen hands over a curve and an accuracy and never sees the epsilon that gets derived. And
bug one is nowhere near that bound: the first curve above runs at nmax = 57.

The fix

Stop when the floats run out. Three lines in the loop:

let r = (scaled_epsilon - 0.5 * (b - a)).max(0.0);
// ... xitp computed exactly as before ...
let xitp = if a < xitp && xitp < b { xitp } else { x1_2 };
if !(a < xitp && xitp < b) {
    return Ok((a, b));
}

Reading them in order.

r is a distance in the paper, so it is never negative there. Unclamped it goes negative as soon as
the bracket falls behind that schedule, and then x1_2 - r.copysign(sigma) moves the point away from
the midpoint instead of toward it, which can put it outside the bracket entirely. Clamping it at
zero is the paper's own behaviour: no slack left, so bisect.

The second line handles a point that rounded onto an endpoint. That can happen with plenty of room
still left in the bracket, since the interpolated point only has to land on a or b by chance.
cubicbez_inv_arclen does exactly that on brackets millions of floats wide. Bisection is this
algorithm's own fallback, so take it rather than give up on a bracket that can still be narrowed.

The third line is the actual exit. If even the midpoint rounds onto an endpoint, there is genuinely
nothing representable left between a and b, every later pass would repeat this one, and this
bracket is the tightest that exists. So return it.

That gives termination without a counter. Each pass either returns or removes at least one
representable value from [a, b], and [a, b] holds a finite number of them.

Build the schedule in floating point. One expression:

-let n1_2 = (((b - a) / epsilon).log2().ceil() - 1.0).max(0.0) as usize;
-let nmax = n0 + n1_2;
-let mut scaled_epsilon = epsilon * (1u64 << nmax) as f64;
+let n1_2 = (((b - a) / epsilon).log2().ceil() - 1.0).max(0.0);
+let nmax = n0 as f64 + n1_2;
+let mut scaled_epsilon = epsilon * nmax.exp2();

nmax comes out of log2 and gets used as 2^nmax. It only ever became an integer so it could
feed a shift. Take the integer out and both overflow sites go with it, since there is no integer
left to overflow.

Nothing here is approximate. f64 holds 2^k exactly for every k below 1024, so exp2(k) and
(1u64 << k) as f64 are the same number for all 64 values the shift could ever have reached.
exp2 needed one line in the define_float_funcs! shim for no_std; libm has exp2 and exp2f.

When nmax now comes out infinite, from a zero or denormal epsilon, scaled_epsilon is NaN,
the clamp on r turns that into a clean zero, and the search runs as plain bisection until the
floats run out. So this commit also drops the 2^-63 lower bound from the doc comment. Nothing
enforces it any more, and epsilon = 0.0 works.

A doc change with no behaviour attached. inv_arclen never said which units accuracy is in.
It is arc length. And a curve of arc length L cannot tell apart arc lengths closer together than
about L * 2^-53, which is roughly 1e-6 on a curve 1e10 long, so asking for 1e-9 there was never
going to work. Both now appear in the doc comment.

Why this is one commit and not two

The two fixes are for independent bugs, and I did write them separately. But the schedule change on
its own would leave things worse than they are today. Without the clamp on r, epsilon = 0.0 gets
you NaN back for a root at 0.9, where today it panics, loudly, in a debug build. Splitting them
puts a commit in the middle that trades a loud failure for a silent one, so I would rather not offer
that as a thing to cherry-pick.

What changes for callers

This does not return the same bits as 0.13.1 on every input, and I would rather say so up front.

I ran 3240 calls through both: 5 brackets, 3 root positions, 6 objective shapes, 9 values of
epsilon, 4 values of n0. Published kurbo returns a root on 1872 of them. This branch gives a
bit-identical answer on 1865 and a different one on 7.

All 7 are the same objective shape, the one whose slope changes by a hundredfold at the root.
Interpolation does badly there, so the bracket falls behind the schedule, r goes negative, and
that is exactly where the clamp changes what happens: published kurbo steps away from the midpoint,
this branch bisects. All 7 answers still land within the epsilon the caller asked for, which is what
solve_itp promises. Some are closer to the true root than before and some are further.

On the 1008 inputs where published kurbo overflows the shift and has no answer at all, this branch
returns one.

That comparison is a program rather than something I am asking you to take on trust. It runs a
transcription of the published loop next to whichever kurbo is linked, and it checks itself against
the crates.io release first so the transcription cannot drift. Happy to share it if you want to run
it yourself.

Tests. 191 kurbo unit tests, 55 doctests and 10 polycool tests all pass. Three new ones cover
the loop ending, a schedule past the width of the shift, and epsilon = 0.0.

The smaller fix I passed on

nmax.min(63) changes one token and needs no shim, so it was the first thing I tried. But on the
1e10 curve it works out scaled_epsilon as 0.9158572095555719, where the right answer,
1.8317144191111439, is a perfectly ordinary f64. The clamp costs one factor of two, so it lands on
exactly half, and the schedule is wrong from the first pass. It also leaves the addition overflow
untouched.

One thing this does not fix

epsilon = f64::NAN makes the loop test false on the very first evaluation, so solve_itp hands
back the midpoint of the bracket without iterating at all. That is 0.5 for a root at 0.9, both
before and after this branch. Same unchecked parameter as the zero case, but narrowing the type of a
public function's argument felt like a bigger call than a bug fix should make on its own.

[mlwilkerson]

This is probably the same -- or overlapping with -- what I reported in #602. And may also overlap with what's in PR #588