From ded35ad2578037cc3d27130d8ae0dd65bbe557cc Mon Sep 17 00:00:00 2001 From: JC Benoist Date: Sat, 21 Mar 2026 13:52:54 +0100 Subject: [PATCH 1/2] Fix bug #1662 --- src/solve.cc | 74 ++++++++++++++++++++++++++++++++++++++-------------- src/tests.cc | 3 +++ 2 files changed, 58 insertions(+), 19 deletions(-) diff --git a/src/solve.cc b/src/solve.cc index c9e3533d5..7f499ae2e 100644 --- a/src/solve.cc +++ b/src/solve.cc @@ -810,6 +810,28 @@ list_p Root::multiple_equation_solver(list_r eqs, list_r names, list_r guesses) } +static array_g jacobi_escape_singular(size_t n, + algebraic_r step, + list_r guesses) +// ---------------------------------------------------------------------------- +// Perturb guesses to escape a degenerate (singular-Jacobian) starting point +// ---------------------------------------------------------------------------- +// Returns the perturbed guesses array, or null on failure. +{ + for (size_t col = 0; col < n; col++) + { + algebraic_g pw = step * decimal::make(col + 1); + if (!pw || !rt.push(+pw)) + return nullptr; + } + array_g perturb = array::from_stack(n, 0); + array_g v = array_p(+guesses); + v = v + perturb; + record(jsolve, "Singular Jacobian, perturbed guesses to %t", +v); + return v; +} + + bool Root::jacobi_solver(list_g &eqs, list_g &vars, list_g &guesses) // ---------------------------------------------------------------------------- // Compute a Jacobian when there is cross-talk between variables @@ -826,9 +848,10 @@ bool Root::jacobi_solver(list_g &eqs, list_g &vars, list_g &guesses) int impr = Settings.SolverImprecision(); algebraic_g eps = algebraic::epsilon(impr); algebraic_g oeps = decimal::make(101,-2); - uint max = Settings.SolverIterations(); - uint iter = 0; - int errs = 0; + uint max = Settings.SolverIterations(); + uint iter = 0; + int errs = 0; + int singular_errs = 0; bool back = false; // Go backwards array_g j, v, d; algebraic_g magnitude, last, forward; @@ -845,20 +868,19 @@ bool Root::jacobi_solver(list_g &eqs, list_g &vars, list_g &guesses) // Set all variables to current value of guesses list::iterator gi = guesses->begin(); + algebraic_g shuffle; + if (errs) + { + // Shuffle slightly around current position (same for all variables) + shuffle = pow(oeps, errs) + eps; + if (!shuffle) + goto error; + } for (object_p varo : *vars) { object_p valo = *gi; - if (errs) - { - if (algebraic_g valg = valo->as_algebraic()) - { - // Shuffle slightly around current position - algebraic_g pw = pow(oeps, errs) + eps; - if (!pw) - goto error; - valo = pw; - } - } + if (errs && valo->as_algebraic()) + valo = shuffle; if (!directory::store_here(varo, valo)) goto error; ++gi; @@ -979,13 +1001,27 @@ bool Root::jacobi_solver(list_g &eqs, list_g &vars, list_g &guesses) j = array::from_stack(n, n, true); v = array::from_stack(n, 0); d = v / j; - record(jsolve, "Jacobian %t values %t delta %t", +j, +v, +d); - v = array_p(+guesses); - v = v - d; + if (!d) + { + if (++singular_errs >= 5) + goto error; + rt.clear_error(); + algebraic_g step = pow(oeps, singular_errs); + if (!step) + goto error; + v = jacobi_escape_singular(n, step, guesses); + last = nullptr; + back = false; + } + else + { + singular_errs = 0; + record(jsolve, "Jacobian %t values %t delta %t", +j, +v, +d); + v = array_p(+guesses); + v = v - d; + } if (!v) goto error; - - // This is the new guesses guesses = +v; } // while (iter < max) diff --git a/src/tests.cc b/src/tests.cc index f26f84387..b002553a5 100644 --- a/src/tests.cc +++ b/src/tests.cc @@ -7744,6 +7744,9 @@ void tests::solver_testing() "{ X Y } { 0 0 } ROOT", ENTER) .expect("{ X=0.5 Y=0.86602 54037 84 }"); + step("Jacobian solver, circle and line with singular initial Jacobian") + .test(CLEAR, "{ 'X^2+Y^2=1' 'X+Y=0' } { X Y } { 0 0 } ROOT", ENTER) + .expect("{ X=-0.70710 67811 87 Y=0.70710 67811 87 }"); step("Solving when the variable is initialized with a constant") .test(CLEAR, DIRECT("m=Ⓒme " "'MSlv(ⒺRelativityMassEnergy;[E];[1 eV])' " From 06adb2154ad4f5cae331a9b4a4c454dfc3936316 Mon Sep 17 00:00:00 2001 From: JC Benoist Date: Sat, 21 Mar 2026 14:28:39 +0100 Subject: [PATCH 2/2] Add another singular jacobian test case --- src/tests.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/tests.cc b/src/tests.cc index 2595a1aea..e90198879 100644 --- a/src/tests.cc +++ b/src/tests.cc @@ -8076,6 +8076,10 @@ void tests::solver_testing() step("Jacobian solver, circle and line with singular initial Jacobian") .test(CLEAR, "{ 'X^2+Y^2=1' 'X+Y=0' } { X Y } { 0 0 } ROOT", ENTER) .expect("{ X=-0.70710 67811 87 Y=0.70710 67811 87 }"); + + step("Jacobian solver, ln equations with singular Jacobian at initial guess") + .test(CLEAR, "{ 'LN(X)+Y=0' 'LN(Y)+X=0' } { X Y } { 1 1 } ROOT", ENTER) + .expect("{ X=0.56714 32904 1 Y=0.56714 32904 1 }"); step("Solving when the variable is initialized with a constant") .test(CLEAR, ("m=Ⓒme " "'MSlv(ⒺRelativityMassEnergy;[E];[1 eV])' "