solver: Add support for multiple equation solving

Add the possibility for the solver to contain a list of equations.
The `NextEq` command rotates the equation list.

Add the `EvalEq` command to evaluate the equation.

Add `EvalEq` and `NextEq` to the solving menu.

Fixes: #1050

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

src/expression.cc

@@ -2130,6 +2130,17 @@ expression_p expression::current_equation(bool error,
         return nullptr;
     }
     id eqty = obj->type();
+    if (eqty == ID_list || eqty == ID_array)
+    {
+        obj = list_p(obj)->at(0);
+        if (!obj)
+        {
+            rt.no_equation_error();
+            return nullptr;
+        }
+        eqty = obj->type();
+    }
+
     if (eqty == ID_equation)
     {
         obj = equation_p(obj)->value();
@@ -2441,6 +2452,24 @@ expression_p expression::as_difference_for_solve() const
 }
 
 
+expression_p expression::left_of_equation() const
+// ----------------------------------------------------------------------------
+//   For the solver, transform A=B into A
+// ----------------------------------------------------------------------------
+{
+    return rewrites(X == Y, X);
+}
+
+
+expression_p expression::right_of_equation() const
+// ----------------------------------------------------------------------------
+//   For the solver, transform A=B into A
+// ----------------------------------------------------------------------------
+{
+    return rewrites(X == Y, Y);
+}
+
+
 expression_p expression::expand() const
 // ----------------------------------------------------------------------------
 //   Run various rewrites to expand terms

src/expression.h

@@ -250,6 +250,8 @@ struct expression : program
     expression_p reorder_terms() const;
     expression_p simplify() const;
     expression_p as_difference_for_solve() const; // Transform A=B into A-B
+    expression_p left_of_equation() const;        // Transform A=B into A
+    expression_p right_of_equation() const;       // Transform A=B into B
     expression_p strip_units(bool keep_constants=false) const;
     object_p     outermost_operator() const;
     size_t       render(renderer &r, bool quoted = false) const

src/ids.tbl

@@ -1024,6 +1024,8 @@ OP(AlgebraVariable,        "ⓧ")        ALIAS(AlgebraVariable,       "VX")
 OP(StoreAlgebraVariable,   "Storeⓧ")   ALIAS(StoreAlgebraVariable,  "StoVX")
 OP(StEq,                   "▶Equation") ALIAS(StEq, "StoreEquation")
 OP(RcEq,                   "Equation▶") ALIAS(RcEq, "RecallEquation")
+OP(NextEq,                 "NextEquation")
+OP(EvalEq,                 "EvaluateEquation")
 
 
 

src/list.cc

@@ -557,7 +557,7 @@ list_p list::append(object_p o) const
 }
 
 
-bool list::expand_without_size() const
+bool list::expand_without_size(size_t *size) const
 // ----------------------------------------------------------------------------
 //   Expand items on the stack, but do not add the size
 // ----------------------------------------------------------------------------
@@ -571,6 +571,8 @@ bool list::expand_without_size() const
             return false;
         }
     }
+    if (size)
+        *size = rt.depth() - depth;
     return true;
 }
 
@@ -818,9 +820,9 @@ HELP_BODY(list)
 //
 // ============================================================================
 
-object::result to_list(uint depth)
+list_p to_list_object(uint depth)
 // ----------------------------------------------------------------------------
-//
+//  Make a list from the stack as an object
 // ----------------------------------------------------------------------------
 {
     scribble scr;
@@ -831,12 +833,27 @@ object::result to_list(uint depth)
             size_t objsz = obj->size();
             byte_p objp = byte_p(obj);
             if (!rt.append(objsz, objp))
-                return object::ERROR;
+                return nullptr;
         }
     }
-    object_g list = list::make(scr.scratch(), scr.growth());
+
-    if (rt.drop(depth) && rt.push(list))
+    if (list_p result = list::make(scr.scratch(), scr.growth()))
-        return object::OK;
+    {
+        rt.drop(depth);
+        return result;
+    }
+    return nullptr;
+}
+
+
+object::result to_list(uint depth)
+// ----------------------------------------------------------------------------
+//  Make a list on the stack
+// ----------------------------------------------------------------------------
+{
+    if (object_g list = to_list_object(depth))
+        if (rt.push(list))
+            return object::OK;
     return object::ERROR;
 }
 

src/list.h

@@ -217,7 +217,7 @@ struct list : text
     }
 
 
-    bool expand_without_size() const;
+    bool expand_without_size(size_t *size = nullptr) const;
     bool expand() const;
     // ------------------------------------------------------------------------
     //   Expand items to the stack, and return number of them
@@ -358,6 +358,7 @@ inline list_g operator*(list_r x, uint y)
 
 
 object::result to_list(uint depth);
+list_p to_list_object(uint depth);
 // ----------------------------------------------------------------------------
 //   Convert `depth` items to a list
 // ----------------------------------------------------------------------------

src/menu.cc

@@ -1197,24 +1197,24 @@ MENU(SolverMenu,
 // ----------------------------------------------------------------------------
 //   The solver menu / application
 // ----------------------------------------------------------------------------
-     "Eq▶",     ID_RcEq,
+     "Eq▶",                     ID_RcEq,
-     "ⓧ",      ID_AlgebraVariable,
+     "ⓧ",                       ID_AlgebraVariable,
-     "Root",    ID_Root,
+     "Root",                    ID_Root,
-     "Solve",   ID_SolvingMenu,
+     "EvalEq",                  ID_EvalEq,
-     "Num",     ID_NumericalSolverMenu,
+     "NxtEq",                   ID_NextEq,
-     "Symb",    ID_SymbolicSolverMenu,
+     "Solve",                   ID_SolvingMenu,
 
-     "▶Eq",     ID_StEq,
+     "▶Eq",                     ID_StEq,
-     "Stoⓧ",    ID_StoreAlgebraVariable,
+     "Stoⓧ",                    ID_StoreAlgebraVariable,
-     "Diff",    ID_DifferentialSolverMenu,
+     "Symb",                    ID_SymbolicSolverMenu,
-     "Poly",    ID_PolynomialSolverMenu,
+     "Diff",                    ID_DifferentialSolverMenu,
-     "Linear",  ID_LinearSolverMenu,
+     "Poly",                    ID_PolynomialSolverMenu,
-     "Multi",   ID_MultiSolverMenu,
+     "Linear",                  ID_LinearSolverMenu,
 
-     "Finance", ID_FinanceSolverMenu,
+     "Multi",                   ID_MultiSolverMenu,
-     "Plot",    ID_PlotMenu,
+     "Finance",                 ID_FinanceSolverMenu,
-     "L.R.",    ID_StatisticsMenu,
+     "Plot",                    ID_PlotMenu,
-     "Eqns",    ID_EquationsMenu,
+     "Eqns",                    ID_EquationsMenu,
      SolverImprecision::label,   ID_SolverImprecision,
      SolverIterations::label,    ID_SolverIterations);
 

src/solve.cc

@@ -31,6 +31,7 @@
 
 #include "algebraic.h"
 #include "arithmetic.h"
+#include "array.h"
 #include "compare.h"
 #include "equations.h"
 #include "expression.h"
@@ -83,8 +84,6 @@ COMMAND_BODY(Root)
         rt.type_error();
         return ERROR;
     }
-    if (eqty == ID_expression)
-        eqobj = expression_p(+eqobj)->as_difference_for_solve();
 
     // Drop input parameters
     rt.drop(3);
@@ -130,6 +129,10 @@ algebraic_p solve(program_g eq, algebraic_g goal, object_g guess)
     object::id  gty = guess->type();
     save<bool>  nodates(unit::nodates, true);
 
+    // Convert A=B+C into A-(B+C)
+    if (eq->type() == object::ID_expression)
+        eq = expression_p(+eq)->as_difference_for_solve();
+
     // Check if low and hight values were given explicitly
     if (gty == object::ID_list || gty == object::ID_array)
     {
@@ -424,9 +427,19 @@ COMMAND_BODY(StEq)
     if (object_p obj = rt.top())
     {
         id objty = obj->type();
+        if (objty == ID_list || objty == ID_array)
+        {
+            for (object_p obj: *list_p(obj))
+            {
+                objty = obj->type();
+                if (objty != ID_expression && objty != ID_polynomial &&
+                    objty != ID_equation)
+                    rt.type_error();
+            }
+        }
         if (objty != ID_expression && objty != ID_polynomial &&
-            objty != ID_equation)
+                objty != ID_equation)
-            rt.type_error();
+                rt.type_error();
         else if (directory::store_here(static_object(ID_Equation), obj))
             if (rt.drop())
                 return OK;
@@ -447,6 +460,80 @@ COMMAND_BODY(RcEq)
 }
 
 
+COMMAND_BODY(NextEq)
+// ----------------------------------------------------------------------------
+//   Cycle equations in the Eq variable if it's a list
+// ----------------------------------------------------------------------------
+{
+    object_p eqname = static_object(ID_Equation);
+    object_p obj = directory::recall_all(eqname, false);
+    if (obj)
+    {
+        id ty = obj->type();
+        if (ty == ID_list || ty == ID_array)
+        {
+            size_t sz = 0;
+            if (list_p(obj)->expand_without_size(&sz))
+            {
+                rt.roll(sz);
+                list_g now = to_list_object(sz);
+                if (directory::store_here(eqname, now))
+                {
+                    ui.menu_refresh(ID_SolvingMenu);
+                    return OK;
+                }
+            }
+        }
+    }
+    return ERROR;
+}
+
+
+COMMAND_BODY(EvalEq)
+// ----------------------------------------------------------------------------
+//   Evaluate the current equation
+// ----------------------------------------------------------------------------
+{
+    if (expression_g expr = expression::current_equation(true))
+    {
+        // We will run programs, do not save stack, etc.
+        settings::PrepareForFunctionEvaluation willEvaluateFunctions;
+
+        expression_g diff = expr->as_difference_for_solve();
+        if (+diff != +expr)
+        {
+            expression_g l = expr->left_of_equation();
+            expression_g r = expr->right_of_equation();
+            algebraic_g lv = l->evaluate();
+            algebraic_g rv = r->evaluate();
+            if (lv && rv)
+            {
+                object_g ltag = tag::make("Left", +lv);
+                object_g rtag = tag::make("Right", +rv);
+                lv = lv - rv;
+                object_g diff = tag::make("Diff", +lv);
+                if (ltag && rtag && diff)
+                {
+                    list_g result = list::make(ID_array, ltag, rtag, diff);
+                    if (result && rt.push(+result))
+                        return OK;
+                }
+            }
+            rt.invalid_equation_error();
+        }
+        else
+        {
+            algebraic_g value = expr->evaluate();
+            if (value)
+                if (tag_p tagged = tag::make("Expr", +value))
+                    if (rt.push(tagged))
+                        return OK;
+        }
+    }
+    return ERROR;
+}
+
+
 MENU_BODY(SolvingMenu)
 // ----------------------------------------------------------------------------
 //   Process the MENU command for SolvingMenu
@@ -455,7 +542,7 @@ MENU_BODY(SolvingMenu)
     expression_p expr   = expression::current_equation(false);
     list_g       vars   = expr ? expr->names() : nullptr;
     size_t       nitems = vars ? vars->items() : 0;
-    items_init(mi, nitems, 3, 1);
+    items_init(mi, nitems+1, 3, 1);
     if (!vars)
         return false;
 
@@ -464,6 +551,7 @@ MENU_BODY(SolvingMenu)
     // First row: Store variables
     mi.plane  = 0;
     mi.planes = 1;
+    menu::items(mi, "EvalEq", ID_EvalEq);
     for (auto name : *vars)
         if (symbol_p sym = name->as<symbol>())
             menu::items(mi, sym, menu::ID_SolvingMenuStore);
@@ -473,6 +561,7 @@ MENU_BODY(SolvingMenu)
     mi.planes = 2;
     mi.skip   = skip;
     mi.index  = mi.plane * ui.NUM_SOFTKEYS;
+    menu::items(mi, "NextEq", ID_NextEq);
     for (auto name : *vars)
         if (symbol_p sym = name->as<symbol>())
             menu::items(mi, sym, menu::ID_SolvingMenuSolve);
@@ -482,6 +571,7 @@ MENU_BODY(SolvingMenu)
     mi.planes = 3;
     mi.skip   = skip;
     mi.index  = mi.plane * ui.NUM_SOFTKEYS;
+    menu::items(mi, "Eq▶", ID_RcEq);
     for (auto name : *vars)
     {
         if (symbol_g sym = name->as<symbol>())
@@ -497,7 +587,7 @@ MENU_BODY(SolvingMenu)
     }
 
     // Add markers
-    for (uint k = 0; k < ui.NUM_SOFTKEYS - (mi.pages > 1); k++)
+    for (uint k = 1; k < ui.NUM_SOFTKEYS - (mi.pages > 1); k++)
     {
         ui.marker(k + 0 * ui.NUM_SOFTKEYS, L'▶', true);
         ui.marker(k + 1 * ui.NUM_SOFTKEYS, L'?', false);
@@ -592,7 +682,7 @@ COMMAND_BODY(SolvingMenuRecall)
     int key = ui.evaluating;
     if (key >= KEY_F1 && key <= KEY_F6)
     {
-        uint index = key - KEY_F1 + 5 * ui.page();
+        uint index = key - KEY_F1 + 5 * ui.page() - 1;
         if (symbol_g sym = expression_variable(index))
             if (object_p value = directory::recall_all(sym, true))
                 if (tag_p tagged = tagged_value(sym, value))
@@ -622,7 +712,7 @@ COMMAND_BODY(SolvingMenuStore)
     int key = ui.evaluating;
     if (key >= KEY_F1 && key <= KEY_F6)
     {
-        uint index = key - KEY_F1 + 5 * ui.page();
+        uint index = key - KEY_F1 + 5 * ui.page() - 1;
         if (algebraic_p entry = expression_variable_or_unit(index))
         {
             if (object_p value = tag::strip(rt.pop()))
@@ -686,7 +776,7 @@ COMMAND_BODY(SolvingMenuSolve)
     int key = ui.evaluating;
     if (key >= KEY_F1 && key <= KEY_F6)
     {
-        uint index = key - KEY_F1 + 5 * ui.page();
+        uint index = key - KEY_F1 + 5 * ui.page() - 1;
         if (symbol_g sym = expression_variable(index))
         {
             object_g value = directory::recall_all(sym, false);

src/solve.h

@@ -40,6 +40,8 @@ COMMAND_DECLARE(Root,3);
 
 COMMAND_DECLARE(StEq, 1);
 COMMAND_DECLARE(RcEq, 0);
+COMMAND_DECLARE(NextEq, 0);
+COMMAND_DECLARE(EvalEq, 0);
 
 
 struct SolvingMenu : menu

state/Demo.48S

@@ -1,3 +1,16 @@
+«
+«
+        0 25 for i
+                i 1_m/s * 1_km/yr convert drop
+        next
+»
+        TEVAL
+»
+
+'UnitBenchmark' STO
+
+
+
 @@ ---------------------------------------------------------------------------
 @@
 @@	 Beep
@@ -706,3 +719,29 @@
         next
 »
 'RandomXYPlot' STO
+
+
+@@ ---------------------------------------------------------------------------
+@@
+@@	 Multiple equations solver
+@@
+@@ ---------------------------------------------------------------------------
+
+{ 'SQ(AA)+SQ(BB)=SQ(CC)'
+
+  'sin(α)/a=sin(β)/b'
+  'sin(β)/b=sin(γ)/c'
+  'sin(α)/b=sin(γ)/c'
+  'α+β+γ=180'
+  's=(a+b+c)/2'
+  'A=sqrt(s*(s-a)*(s-b)*(s-c))'
+  'sq(c)=sq(a)+sq(b)-2*a*b*cos(γ)'
+  'sq(b)=sq(a)+sq(c)-2*a*c*cos(β)'
+  'sq(a)=sq(b)+sq(c)-2*b*c*cos(α)'
+} STEQ
+
+10_m 'AA' STO
+85_ft 'BB' STO
+1_cm 'CC' STO
+
+SolvingMenu