Solver fails to find solution

[c3d]

{ 'X^2+Y^2=1' 'X+Y=0' } { X Y } { 0 0 } ROOT

The solver fails stating it cannot find a solution. It works with X-Y=0 as the second equation.

[jcbenoist]

I am working on it (well Claude also, to be honest), and this will teach me the algebra section of db48x I do not know for now.

[jcbenoist]

Here's a summary of Claude investigations. I'll try to fully understand it before publishing (otherwise it's too easy and uninteresting), and let first comments arrives

Root cause: At initial guess (0,0), the Jacobian for {X²+Y²=1, X+Y=0} is singular because ∂(X+Y)/∂X = ∂(X+Y)/∂Y = 1. The matrix division v/j called j->invert() which detected zero determinant, set rt.zero_divide_error(), and returned null. The runtime error state then blocked all subsequent operations.

Fix in [src/solve.cc]:

Added int singular_errs = 0 counter
After d = v / j, when d is null: clear the runtime error (rt.clear_error()), increment singular_errs, push column-dependent perturbations [step, 2*step, ...] where step = oeps^singular_errs, add them to guesses, reset last = nullptr and back = false (to avoid using the now-invalid d in the next iteration's direction check), then continue
After 5 consecutive singular Jacobians → give up

New test in [src/tests.cc] — line 7747: { 'X^2+Y^2=1' 'X+Y=0' } { X Y } { 0 0 } ROOT → { X=-0.70710⁠67811⁠87 Y=0.70710⁠67811⁠87 } (converges to the -X, +Y solution branch from the perturbation at [1.01, 2.02]).