integrate: Initial implementation of numerical integration

Implement a really basic algorithm for integration evaluation.
This divides the integration surface by two every time.
This appears much less efficient than the integration done by HP calculators.

Fixes: #42

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

Makefile

@@ -213,6 +213,7 @@ CXX_SOURCES +=				\
 	src/grob.cc			\
 	src/plot.cc			\
 	src/solve.cc			\
+	src/integrate.cc		\
 	fonts/HelpFont.cc		\
 	fonts/EditorFont.cc		\
 	fonts/StackFont.cc

fonts/EditorFont.cc

@@ -7,7 +7,7 @@ extern const unsigned char EditorFont_sparse_font_data[];
 const unsigned char EditorFont_sparse_font_data[73676] FONT_QSPI =
 {
 
-    0xD7, 0x02, 0xC7, 0xBF, 0x04, 0x36, 0x00, 0x01, 0x00, 0x2C, 0x01, 0x01, 0x00, 0x00, 0x0D, 0x01,
+    0xD8, 0x02, 0xC7, 0xBF, 0x04, 0x36, 0x00, 0x01, 0x00, 0x2C, 0x01, 0x01, 0x00, 0x00, 0x0D, 0x01,
     0x00, 0x2C, 0x01, 0x01, 0x08, 0x00, 0x20, 0x5F, 0x00, 0x2C, 0x01, 0x01, 0x08, 0x00, 0x02, 0x0B,
     0x05, 0x21, 0x09, 0xEF, 0xBD, 0xF7, 0xDE, 0x7B, 0xEF, 0xBD, 0xF7, 0xDE, 0x7B, 0xEF, 0xBD, 0xF7,
     0x1E, 0x00, 0x00, 0xB8, 0xFF, 0xFF, 0xFF, 0x0E, 0x03, 0x0B, 0x0B, 0x0F, 0x11, 0x8F, 0x7F, 0xFC,

fonts/HelpFont.cc

@@ -7,7 +7,7 @@ extern const unsigned char HelpFont_sparse_font_data[];
 const unsigned char HelpFont_sparse_font_data[15010] FONT_QSPI =
 {
 
-    0xD7, 0x02, 0x9E, 0x75, 0x14, 0x00, 0x01, 0x00, 0x11, 0x01, 0x01, 0x00, 0x00, 0x0D, 0x01, 0x00,
+    0xD8, 0x02, 0x9E, 0x75, 0x14, 0x00, 0x01, 0x00, 0x11, 0x01, 0x01, 0x00, 0x00, 0x0D, 0x01, 0x00,
     0x11, 0x01, 0x01, 0x03, 0x00, 0x20, 0x5F, 0x00, 0x11, 0x01, 0x01, 0x03, 0x00, 0x01, 0x04, 0x02,
     0x0D, 0x04, 0xFF, 0xFF, 0xF3, 0x03, 0x01, 0x04, 0x04, 0x05, 0x06, 0xFF, 0xFF, 0x0F, 0x01, 0x04,
     0x08, 0x0D, 0x0A, 0x12, 0x12, 0x14, 0x7F, 0x7F, 0x24, 0x24, 0x24, 0xFE, 0xFE, 0x28, 0x48, 0x48,

fonts/StackFont.cc

@@ -7,7 +7,7 @@ extern const unsigned char StackFont_sparse_font_data[];
 const unsigned char StackFont_sparse_font_data[36409] FONT_QSPI =
 {
 
-    0xD7, 0x02, 0xB4, 0x9C, 0x02, 0x24, 0x00, 0x01, 0x00, 0x1C, 0x01, 0x01, 0x00, 0x00, 0x0D, 0x01,
+    0xD8, 0x02, 0xB4, 0x9C, 0x02, 0x24, 0x00, 0x01, 0x00, 0x1C, 0x01, 0x01, 0x00, 0x00, 0x0D, 0x01,
     0x00, 0x1C, 0x01, 0x01, 0x06, 0x00, 0x20, 0x5F, 0x00, 0x1C, 0x01, 0x01, 0x06, 0x00, 0x02, 0x05,
     0x03, 0x17, 0x06, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x1F, 0xC0, 0xFF, 0x1F, 0x02, 0x05, 0x08, 0x0A,
     0x0C, 0xE7, 0xE7, 0xE7, 0xE7, 0xE7, 0xE7, 0xE7, 0xE7, 0xE7, 0xE7, 0x01, 0x05, 0x0E, 0x17, 0x10,

sim/simulator.pro

@@ -79,6 +79,7 @@ SOURCES +=                                      \
 	../src/grob.cc                          \
 	../src/plot.cc	                        \
 	../src/solve.cc	                        \
+	../src/integrate.cc                     \
 	../src/tests.cc
 
 HEADERS +=                                      \

src/errors.tbl

@@ -65,8 +65,7 @@ ERROR(number_too_big,           "Number is too big")
 ERROR(too_many_rewrites,        "Too many rewrites")
 ERROR(invalid_ppar,             "Invalid PlotParameters")
 ERROR(invalid_plot_type,        "Invalid plot type")
-ERROR(invalid_plot_function,    "Invalid plot function")
+ERROR(invalid_function,         "Invalid function")
-ERROR(invalid_solve_function,   "Invalid solve function")
 ERROR(constant_value,           "Constant?")
 ERROR(bad_guess,                "Bad guess?")
 ERROR(no_solution,              "No solution?")

src/ids.tbl

@@ -385,9 +385,9 @@ CMD(Type)
 CMD(TypeName)
 CMD(Version)
 
-// Numerical solver
+// High-level applications
 CMD(Root)
-
+NAMED(Integrate, "∫")
 
 // ============================================================================
 //

src/integrate.cc

@@ -0,0 +1,203 @@
+// ****************************************************************************
+//  integrate.cc                                                 DB48X project
+// ****************************************************************************
+//
+//   File Description:
+//
+//
+//
+//
+//
+//
+//
+//
+//
+//
+// ****************************************************************************
+//   (C) 2023 Christophe de Dinechin <christophe@dinechin.org>
+//   This software is licensed under the terms outlined in LICENSE.txt
+// ****************************************************************************
+//   This file is part of DB48X.
+//
+//   DB48X is free software: you can redistribute it and/or modify
+//   it under the terms outlined in the LICENSE.txt file
+//
+//   DB48X is distributed in the hope that it will be useful,
+//   but WITHOUT ANY WARRANTY; without even the implied warranty of
+//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// ****************************************************************************
+
+#include "integrate.h"
+
+#include "algebraic.h"
+#include "arithmetic.h"
+#include "compare.h"
+#include "equation.h"
+#include "functions.h"
+#include "integer.h"
+#include "recorder.h"
+#include "settings.h"
+#include "symbol.h"
+#include "tag.h"
+
+RECORDER(integrate,         16, "Numerical integration");
+RECORDER(integrate_error,   16, "Numerical integrationsol");
+
+
+COMMAND_BODY(Integrate)
+// ----------------------------------------------------------------------------
+//   Numerical integration
+// ----------------------------------------------------------------------------
+{
+    if (!rt.args(4))
+        return ERROR;
+
+    object_g variable = rt.stack(0);
+    object_g eq = rt.stack(1);
+    object_g high = rt.stack(2);
+    object_g low = rt.stack(3);
+    if (!eq || !variable || !high || !low)
+        return ERROR;
+
+    record(integrate,
+           "Integrating %t for variable %t in range %t-%t",
+           eq.Safe(), variable.Safe(), low.Safe(), high.Safe());
+
+    // Check that we have a variable name on stack level 1 and
+    // a proram or equation on level 2
+    symbol_g name = variable->as_quoted<symbol>();
+    id eqty = eq->type();
+    if (eqty != ID_program && eqty != ID_equation)
+        name = nullptr;
+    if (!name || !low->is_algebraic() || !high->is_algebraic())
+    {
+        rt.type_error();
+        return ERROR;
+    }
+
+    // Drop input parameters
+    rt.drop(4);
+
+    // Actual integration
+    algebraic_g intg = integrate(eq, name,
+                                 algebraic_p(low.Safe()),
+                                 algebraic_p(high.Safe()));
+    if (intg.Safe() &&rt.push(intg.Safe()))
+        return OK;
+
+    return ERROR;
+}
+
+
+algebraic_p integrate(object_g     eq,
+                      symbol_g     name,
+                      algebraic_g lx,
+                      algebraic_g hx)
+// ----------------------------------------------------------------------------
+//   The core of the integration function
+// ----------------------------------------------------------------------------
+{
+    // Check if the guess is an algebraic or if we need to extract one
+    algebraic_g x, dx, dx2;
+    algebraic_g y, dy, sy, sy2;
+    algebraic_g two = integer::make(2);
+    record(integrate, "Initial range %t-%t", lx.Safe(), hx.Safe());
+
+    // Set independent variable
+    save<symbol_g *> iref(equation::independent, &name);
+    save<object_g *> ival(equation::independent_value, (object_g *) &x);
+    int              prec = -Settings.integprec;
+    algebraic_g      eps = rt.make<decimal128>(object::ID_decimal128,
+                                               prec, true);
+
+    // Initial integration step and sum
+    dx = hx - lx;
+    sy = integer::make(0);
+    if (!dx || !sy)
+        return nullptr;
+
+    // Loop for a maximum number of conversion iterations
+    uint max = Settings.maxinteg;
+    uint iter = 1;
+    for (uint d = 0; iter <= max && !program::interrupted(); d++)
+    {
+        sy2 = sy;
+        dx2 = dx / two;
+        x   = lx + dx2;
+        sy  = integer::make(0);
+        if (!x || !sy)
+            return nullptr;
+
+        for (uint i = 0; i < iter; i++)
+        {
+            // If we are starting to use really big numbers, approximate
+            if (x->is_big())
+                if (!algebraic::to_decimal(x))
+                    return nullptr;
+
+            // Evaluate equation
+            size_t depth = rt.depth();
+            if (!rt.push(x.Safe()))
+                return nullptr;
+            record(integrate, "[%u:%u] x=%t", d, i, x.Safe());
+
+            object::result err    = eq->execute();
+            size_t         dnow   = rt.depth();
+            if (dnow != depth + 1 && dnow != depth + 2)
+            {
+                record(integrate_error, "Depth moved from %u to %u", depth, dnow);
+                rt.invalid_function_error();
+                return nullptr;
+            }
+
+            if (err != object::OK)
+            {
+                // Error on last function evaluation, try again
+                record(integrate_error, "Got error %+s", rt.error());
+                return nullptr;
+            }
+
+            y = algebraic_p(rt.pop());
+            if (dnow == depth + 2)
+                rt.drop();
+            record(integrate, "[%u:%u] x=%t y=%t", d, i, x.Safe(), y.Safe());
+
+            if (!y)
+                return nullptr;
+            if (!y->is_algebraic())
+            {
+                rt.invalid_function_error();
+                return nullptr;
+            }
+
+            sy = sy + y;
+            x = x + dx;
+            record(integrate, "[%u:%u] sy=%t", d, i, sy.Safe());
+            if (!sy || !x)
+                return nullptr;
+
+            if (sy->is_big())
+                if (!algebraic::to_decimal(sy))
+                    return nullptr;
+        }
+        sy = sy * dx;
+        record(integrate, "[%u] Sum sy=%t sy2=%t dx=%t",
+               d, sy.Safe(), sy2.Safe(), dx.Safe());
+        if (!sy)
+            return nullptr;
+        if (smaller_magnitude(sy - sy2, eps * sy2))
+        {
+            sy = (sy + sy2) / two;
+            break;
+        }
+
+        dx = dx2;
+        sy = sy + sy2 / two;
+        if (!sy)
+            return nullptr;
+
+        iter <<= 1;
+    }
+
+    return sy;
+}

src/integrate.h

@@ -0,0 +1,43 @@
+#ifndef INTEGRATE_H
+#define INTEGRATE_H
+// ****************************************************************************
+//  integrate.h                                                       DB48X project
+// ****************************************************************************
+//
+//   File Description:
+//
+//
+//
+//
+//
+//
+//
+//
+//
+//
+// ****************************************************************************
+//   (C) 2023 Christophe de Dinechin <christophe@dinechin.org>
+//   This software is licensed under the terms outlined in LICENSE.txt
+// ****************************************************************************
+//   This file is part of DB48X.
+//
+//   DB48X is free software: you can redistribute it and/or modify
+//   it under the terms outlined in the LICENSE.txt file
+//
+//   DB48X is distributed in the hope that it will be useful,
+//   but WITHOUT ANY WARRANTY; without even the implied warranty of
+//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// ****************************************************************************
+
+#include "algebraic.h"
+#include "command.h"
+#include "symbol.h"
+
+algebraic_p integrate(object_g    eq,
+                      symbol_g    name,
+                      algebraic_g low,
+                      algebraic_g high);
+
+COMMAND_DECLARE(Integrate);
+
+#endif // INTEGRATE_H

src/menu.cc

@@ -875,7 +875,7 @@ MENU(IntegrationMenu,
 // ----------------------------------------------------------------------------
 //   Symbolic and numerical integration
 // ----------------------------------------------------------------------------
-     "Num",     ID_Unimplemented,
+     "Num",     ID_Integrate,
      "Symb",    ID_Unimplemented,
      "Prim",    ID_Unimplemented,
 

src/object.cc

@@ -48,6 +48,7 @@
 #include "graphics.h"
 #include "grob.h"
 #include "integer.h"
+#include "integrate.h"
 #include "list.h"
 #include "locals.h"
 #include "logical.h"

src/plot.cc

@@ -188,7 +188,7 @@ object::result draw_plot(object::id            kind,
         else
         {
             if (!rt.error())
-                rt.invalid_plot_function_error();
+                rt.invalid_function_error();
             if (rt.depth() > depth)
                 rt.drop(rt.depth() - depth);
             Screen.text(0, 0, rt.error(), ErrorFont,

src/settings.h

@@ -81,7 +81,9 @@ struct settings
           wordsize(64),
           maxbignum(1024),
           maxsolve(1024),
+          maxinteg(1024),
           solveprec(24),
+          integprec(24),
           maxrewrites(100),
           fraciter(10),
           fracprec(12),
@@ -189,7 +191,9 @@ public:
     size_t   wordsize;          // Wordsize for binary numbers (in bits)
     size_t   maxbignum;         // Maximum size for a bignum (in bits)
     size_t   maxsolve;          // Maximum number of iterations for solver
+    size_t   maxinteg;          // Maximum number of iterations for integration
     uint16_t solveprec;         // Precision of solver in digits
+    uint16_t integprec;         // Precision of integration in digits
     uint16_t maxrewrites;       // Maximum number of rewrites
     uint16_t fraciter;          // Number of iterations for ->Q
     uint16_t fracprec;          // Number of digits for ->Q

src/solve.cc

@@ -144,7 +144,7 @@ algebraic_p solve(object_g eq, symbol_g name, object_g guess)
         if (dnow != depth + 1 && dnow != depth + 2)
         {
             record(solve_error, "Depth moved from %u to %u", depth, dnow);
-            rt.invalid_solve_function_error();
+            rt.invalid_function_error();
             return nullptr;
         }
         if (err != object::OK)
@@ -165,9 +165,11 @@ algebraic_p solve(object_g eq, symbol_g name, object_g guess)
             if (dnow == depth + 2)
                 rt.drop();
             record(solve, "[%u] x=%t y=%t", i, x.Safe(), y.Safe());
-            if (!y || !y->is_algebraic())
+            if (!y)
+                return nullptr;
+            if (!y->is_algebraic())
             {
-                rt.invalid_solve_function_error();
+                rt.invalid_function_error();
                 return nullptr;
             }
             if (y->is_zero() || smaller_magnitude(y, eps))
@@ -278,7 +280,7 @@ algebraic_p solve(object_g eq, symbol_g name, object_g guess)
             // Check if there are unresolved symbols
             if (x->is_strictly_symbolic())
             {
-                rt.invalid_solve_function_error();
+                rt.invalid_function_error();
                 return x;
             }
 
@@ -287,7 +289,7 @@ algebraic_p solve(object_g eq, symbol_g name, object_g guess)
             {
                 if (!algebraic::to_decimal(x))
                 {
-                    rt.invalid_solve_function_error();
+                    rt.invalid_function_error();
                     return x;
                 }
             }
@@ -318,7 +320,7 @@ algebraic_p solve(object_g eq, symbol_g name, object_g guess)
            x.Safe(), y.Safe(), lx.Safe(), ly.Safe());
 
     if (!is_valid)
-        rt.invalid_solve_function_error();
+        rt.invalid_function_error();
     else if (is_constant)
         rt.constant_value_error();
     else