solver: Adjust detection of "epsilon" for large values

Compute the "epsilon" for X and Y axis as follows:

- For the X axis, use the initial range of X
- For the Y axis, use the first computed value for the function

This fixes the original problem in #1179, where the equation was in
the form `f(x)=0`, and we had a heuristic in that case that would
deduce an incorrect epsilon from the right-hand-side value `0`.

That heuristic was implemented to solve equations with very large
values, such as found in physics when dealing for example with
Avogadro number. A simple example would be something like
`1E45*sin(x)=0.5E45`. But the heuristic was depending on being able to
isolate the left and right hand sides. It would fail for example if
the expresion was initially written as `1E45*sin(x)-0.5E45`.

Instead, the new heuristic computes the epsilon for the Y axis from
the first run of the equation. In that case, if the original epsilon
is `1E-18` (which is the default for a `Precision` of 24 digits and a
`SolverImprecision` value of 6), then the target precision will be in
the order of `1E45*1E-18`.

Fixes: #1179

Signed-off-by: Christophe de Dinechin <christophe@dinechin.org>

src/solve.cc

@@ -127,19 +127,10 @@ algebraic_p solve(program_g eq, algebraic_g goal, object_g guess)
     save<bool>  nodates(unit::nodates, true);
 
     // Convert A=B+C into A-(B+C)
-    program_g left, right;
-    bool iseq = false;
     if (eq->type() == object::ID_expression)
-    {
-        expression_g diff = expression_p(+eq)->as_difference_for_solve();
-        if (+diff != +eq)
-        {
-            left = expression_p(+eq)->left_of_equation();
-            right = expression_p(+eq)->right_of_equation();
-            eq = +diff;
-            iseq = left && right;
-        }
-    }
+        if (expression_g diff = expression_p(+eq)->as_difference_for_solve())
+            if (+diff != +eq)
+                eq = +diff;
 
     // Check if low and hight values were given explicitly
     if (gty == object::ID_list || gty == object::ID_array)
@@ -232,7 +223,8 @@ algebraic_p solve(program_g eq, algebraic_g goal, object_g guess)
     }
     save<symbol_g *> iref(expression::independent, &name);
     int              prec = Settings.Precision() - Settings.SolverImprecision();
-    algebraic_g      eps  = decimal::make(1, prec <= 0 ? -1 : -prec);
+    algebraic_g      yeps  = decimal::make(1, prec <= 0 ? -1 : -prec);
+    algebraic_g      xeps  = (lx + hx) * yeps;
     bool             is_constant = true;
     bool             is_valid    = false;
     uint             max         = Settings.SolverIterations();
@@ -246,23 +238,21 @@ algebraic_p solve(program_g eq, algebraic_g goal, object_g guess)
             return x;
 
         // Evaluate equation
-        if (iseq)
+        y = algebraic::evaluate_function(eq, x);
+
+        // If the function evaluates as 10^23 and eps=10^-18, use 10^(23-18)
+        if (!i && y && !y->is_zero())
         {
-            dx = algebraic::evaluate_function(left, x);
-            dy = algebraic::evaluate_function(right, x);
-            y = dx - dy;
-            dx = dy * eps;
-            if (dx)
-                if (unit_p ru = dx->as<unit>())
-                    dx = ru->value();
-        }
-        else
-        {
-            y = algebraic::evaluate_function(eq, x);
-            dx = eps;
+            if (algebraic_g neps = abs::run(y) * yeps)
+            {
+                if (unit_p ru = neps->as<unit>())
+                    neps = ru->value();
+                if (smaller_magnitude(yeps, neps))
+                    yeps = neps;
+            }
         }
         record(solve, "[%u] x=%t (%t to %t)  y=%t (%t to %t) err=%t",
-               i, +x, +lx, +hx, +y, +ly, +hy, +dx);
+               i, +x, +lx, +hx, +y, +ly, +hy, +yeps);
         if (!y)
         {
             // Error on last function evaluation, try again
@@ -283,7 +273,7 @@ algebraic_p solve(program_g eq, algebraic_g goal, object_g guess)
                 dy = yu->value();
             else
                 dy = y;
-            if (dy->is_zero() || smaller_magnitude(dy, dx))
+            if (dy->is_zero() || smaller_magnitude(dy, yeps))
             {
                 record(solve, "[%u] Solution=%t value=%t", i, +x, +y);
                 return x;
@@ -346,10 +336,10 @@ algebraic_p solve(program_g eq, algebraic_g goal, object_g guess)
                 return nullptr;
             dy = hx + lx;
             if (!dy || dy->is_zero(false))
-                dy = eps;
+                dy = yeps;
             else
-                dy = dy * eps;
-            if (dx->is_zero(false) || smaller_magnitude(dx, dy))
+                dy = dy * yeps;
+            if (dx->is_zero(false) || smaller_magnitude(dx, xeps))
             {
                 x = lx;
                 record(solve, "[%u] Minimum=%t value=%t", i, +x, +y);
@@ -412,7 +402,7 @@ algebraic_p solve(program_g eq, algebraic_g goal, object_g guess)
                                      object::ID_Deg);
                 else
                     dx = integer::make(0x1081 * s * i);
-                dx = dx * eps;
+                dx = dx * yeps;
                 if (x && x->is_zero())
                     x = dx;
                 else

src/tests.cc

@@ -172,7 +172,7 @@ void tests::run(uint onlyCurrent)
     {
         here().begin("Current");
         if (onlyCurrent & 1)
-            editor_operations();
+            solver_testing();
         if (onlyCurrent & 2)
             demo_ui();
         if (onlyCurrent & 4)
@@ -5735,6 +5735,13 @@ void tests::solver_testing()
         .test(LSHIFT, A, LSHIFT, A)
         .expect("0.5 m");
 
+    step("Solving with large values (#1179")
+        .test(CLEAR, "DEG '1E45*sin(x)-0.5E45' 'x' 2 ROOT", ENTER)
+        .expect("x:30.");
+    step("Solving equation containing a zero side (#1179")
+        .test(CLEAR, "'-3*expm1(-x)-x=0' 'x' 2 ROOT", ENTER)
+        .expect("x:2.82143 93721 2");
+
     step("Exit: Clear variables")
         .test(CLEAR, "UPDIR 'SLVTST' PURGE", ENTER);
 }
@@ -5807,7 +5814,7 @@ void tests::eqnlib_columns_and_beams()
         .test(CLEAR, LSHIFT, F1, LSHIFT, F4)
         .expect("r:4.1148 cm")
         .test(NOSHIFT, F1)
-        .expect("'4.1148 cm=411.48 mm↑2/cm+1.212⁳⁻¹⁸ mm↑2/cm'");
+        .expect("'4.1148 cm=411.48 mm↑2/cm+6.12⁳⁻¹⁹ mm↑2/cm'");
 
     step("Solving Eccentric Columns")
         .test(CLEAR, RSHIFT, F, F2, RSHIFT, F2)
@@ -10865,10 +10872,8 @@ tests &tests::editing(size_t length, uint extrawait)
 {
     nokeys(extrawait);
     return check(rt.editing() == length,
-                 "Expected editing length to be ",
-                 length,
-                 " got ",
-                 rt.editing());
+                 "Expected editing length to be ", length,
+                 " got ", rt.editing());
 }
 
 
@@ -10888,12 +10893,8 @@ tests &tests::editor(cstring text, uint extrawait)
     {
         if (rt.error())
         {
-            explain("Expected editor [",
-                    text,
-                    "], "
-                    "got error [",
-                    rt.error(),
-                    "] instead");
+            explain("Expected editor [", text, "], "
+                    "got error [", rt.error(), "] instead");
             return fail();
         }
 
@@ -10906,30 +10907,16 @@ tests &tests::editor(cstring text, uint extrawait)
     }
 
     if (!ed)
-        explain("Expected editor to contain [",
-                text,
-                "], "
+        explain("Expected editor to contain [", text, "], "
                 "but it's empty");
     if (sz != strlen(text))
-        explain("Expected ",
-                strlen(text),
-                " characters in editor"
-                " [",
-                text,
-                "], "
-                "but got ",
-                sz,
-                " characters "
-                " [",
-                std::string(cstring(ed), sz),
-                "]");
+        explain("Expected ", strlen(text), " characters in editor"
+                " [", text, "], "
+                "but got ", sz, " characters "
+                " [", std::string(cstring(ed), sz), "]");
     if (memcmp(ed, text, sz))
-        explain("Expected editor to contain [",
-                text,
-                "], "
-                "but it contains [",
-                std::string(cstring(ed), sz),
-                "]");
+        explain("Expected editor to contain [", text, "], "
+                "but it contains [", std::string(cstring(ed), sz), "]");
 
     fail();
     return *this;