From 6f3c07e759d9770edf8a12148fd2f85330dd05b2 Mon Sep 17 00:00:00 2001 From: Christophe de Dinechin Date: Thu, 12 Oct 2023 22:59:20 +0200 Subject: [PATCH] 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 --- Makefile | 1 + fonts/EditorFont.cc | 2 +- fonts/HelpFont.cc | 2 +- fonts/StackFont.cc | 2 +- sim/simulator.pro | 1 + src/errors.tbl | 3 +- src/ids.tbl | 4 +- src/integrate.cc | 203 ++++++++++++++++++++++++++++++++++++++++++++ src/integrate.h | 43 ++++++++++ src/menu.cc | 2 +- src/object.cc | 1 + src/plot.cc | 2 +- src/settings.h | 4 + src/solve.cc | 14 +-- 14 files changed, 269 insertions(+), 15 deletions(-) create mode 100644 src/integrate.cc create mode 100644 src/integrate.h diff --git a/Makefile b/Makefile index d1e7a974..50d748b7 100644 --- a/Makefile +++ b/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 diff --git a/fonts/EditorFont.cc b/fonts/EditorFont.cc index 298126d7..fa146252 100644 --- a/fonts/EditorFont.cc +++ b/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, diff --git a/fonts/HelpFont.cc b/fonts/HelpFont.cc index d839ca54..a1e74ec7 100644 --- a/fonts/HelpFont.cc +++ b/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, diff --git a/fonts/StackFont.cc b/fonts/StackFont.cc index c0d29eee..ecfe6ce8 100644 --- a/fonts/StackFont.cc +++ b/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, diff --git a/sim/simulator.pro b/sim/simulator.pro index 4081f140..d1ecaada 100644 --- a/sim/simulator.pro +++ b/sim/simulator.pro @@ -79,6 +79,7 @@ SOURCES += \ ../src/grob.cc \ ../src/plot.cc \ ../src/solve.cc \ + ../src/integrate.cc \ ../src/tests.cc HEADERS += \ diff --git a/src/errors.tbl b/src/errors.tbl index d4e7ffd0..78af1cb9 100644 --- a/src/errors.tbl +++ b/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_solve_function, "Invalid solve function") +ERROR(invalid_function, "Invalid function") ERROR(constant_value, "Constant?") ERROR(bad_guess, "Bad guess?") ERROR(no_solution, "No solution?") diff --git a/src/ids.tbl b/src/ids.tbl index a41794c0..b2ca548f 100644 --- a/src/ids.tbl +++ b/src/ids.tbl @@ -385,9 +385,9 @@ CMD(Type) CMD(TypeName) CMD(Version) -// Numerical solver +// High-level applications CMD(Root) - +NAMED(Integrate, "∫") // ============================================================================ // diff --git a/src/integrate.cc b/src/integrate.cc new file mode 100644 index 00000000..deb521a9 --- /dev/null +++ b/src/integrate.cc @@ -0,0 +1,203 @@ +// **************************************************************************** +// integrate.cc DB48X project +// **************************************************************************** +// +// File Description: +// +// +// +// +// +// +// +// +// +// +// **************************************************************************** +// (C) 2023 Christophe de Dinechin +// 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(); + 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 iref(equation::independent, &name); + save ival(equation::independent_value, (object_g *) &x); + int prec = -Settings.integprec; + algebraic_g eps = rt.make(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; +} diff --git a/src/integrate.h b/src/integrate.h new file mode 100644 index 00000000..978597e7 --- /dev/null +++ b/src/integrate.h @@ -0,0 +1,43 @@ +#ifndef INTEGRATE_H +#define INTEGRATE_H +// **************************************************************************** +// integrate.h DB48X project +// **************************************************************************** +// +// File Description: +// +// +// +// +// +// +// +// +// +// +// **************************************************************************** +// (C) 2023 Christophe de Dinechin +// 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 diff --git a/src/menu.cc b/src/menu.cc index 1cd539b4..024df2f5 100644 --- a/src/menu.cc +++ b/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, diff --git a/src/object.cc b/src/object.cc index d0421255..294f0dc2 100644 --- a/src/object.cc +++ b/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" diff --git a/src/plot.cc b/src/plot.cc index b5f6c9e6..8d98fb6f 100644 --- a/src/plot.cc +++ b/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, diff --git a/src/settings.h b/src/settings.h index 924520ad..148803cb 100644 --- a/src/settings.h +++ b/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 diff --git a/src/solve.cc b/src/solve.cc index 7157cb4c..3f0b10a7 100644 --- a/src/solve.cc +++ b/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