lvgl_micropython.diff

diff -urN lvgl_micropython.O/ext_mod/micropython.cmake lvgl_micropython/ext_mod/micropython.cmake
--- lvgl_micropython.O/ext_mod/micropython.cmake	2024-11-25 14:57:11.128573936 +0100
+++ lvgl_micropython/ext_mod/micropython.cmake	2024-11-25 13:49:48.934368008 +0100
@@ -1,4 +1,5 @@
 include(${CMAKE_CURRENT_LIST_DIR}/lcd_bus/micropython.cmake)
 include(${CMAKE_CURRENT_LIST_DIR}/lvgl/micropython.cmake)
 include(${CMAKE_CURRENT_LIST_DIR}/lcd_utils/micropython.cmake)
+include(${CMAKE_CURRENT_LIST_DIR}/ulab/code/micropython.cmake)
 
diff -urN lvgl_micropython.O/ext_mod/ulab/code/micropython.mk lvgl_micropython/ext_mod/ulab/code/micropython.mk
--- lvgl_micropython.O/ext_mod/ulab/code/micropython.mk	2024-11-25 14:59:24.496963159 +0100
+++ lvgl_micropython/ext_mod/ulab/code/micropython.mk	2024-11-25 14:21:51.078125041 +0100
@@ -2,6 +2,7 @@
 USERMODULES_DIR := $(USERMOD_DIR)
 
 # Add all C files to SRC_USERMOD.
+SRC_USERMOD += $(USERMODULES_DIR)/scipy/integrate/integrate.c
 SRC_USERMOD += $(USERMODULES_DIR)/scipy/linalg/linalg.c
 SRC_USERMOD += $(USERMODULES_DIR)/scipy/optimize/optimize.c
 SRC_USERMOD += $(USERMODULES_DIR)/scipy/signal/signal.c
diff -urN lvgl_micropython.O/ext_mod/ulab/code/scipy/integrate/integrate.c lvgl_micropython/ext_mod/ulab/code/scipy/integrate/integrate.c
--- lvgl_micropython.O/ext_mod/ulab/code/scipy/integrate/integrate.c	1970-01-01 01:00:00.000000000 +0100
+++ lvgl_micropython/ext_mod/ulab/code/scipy/integrate/integrate.c	2024-11-25 14:47:40.946307334 +0100
@@ -0,0 +1,635 @@
+/*
+ * This file is not part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2024 Harald Milz <hm@seneca.muc.de>
+ *
+ * References:
+ * - Dr. Robert van Engelen, Improving the mp_float_t Exponential Quadrature Tanh-Sinh, Sinh-Sinh and Exp-Sinh Formulas,
+ *   2021, https://www.genivia.com/qthsh.html 
+ * - Borwein, Bailey & Girgensohn, "Experimentation in Mathematics - Computational Paths to Discovery", A K Peters,
+ *   2003, pages 312-313
+ * - Joren Vanherck, Bart Sorée, Wim Magnus, Tanh-sinh quadrature for single and multiple integration using 
+ *   floating-point arithmetic, 2020, https://arxiv.org/abs/2007.15057
+ * - Tanh-Sinh quadrature, Wikipedia, https://en.wikipedia.org/wiki/Tanh-sinh_quadrature
+ * - Romberg's method, Wikipedia, https://en.wikipedia.org/wiki/Romberg%27s_method
+ * - Adaptive Simpson's method, Wikipedia, https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method
+ * - Gauss–Kronrod quadrature formula, Wikipedia, https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula
+ *  
+ * This module provides four integration methods, and thus deviates from scipy.integrate a bit. 
+ * As for the pros and cons of the different methods please consult the literature above. 
+ * The code was ported to Micropython from Dr. Engelen's paper and used with his written kind permission
+ * - quad    - Tanh-Sinh, Sinh-Sinh and Exp-Sinh quadrature
+ * - romberg - Romberg quadrature
+ * - simpson - Adaptive Simpson quadrature
+ * - quadgk  - Adaptive Gauss-Kronrod (G10,K21) quadrature
+ */
+
+#include <math.h>
+#include "py/obj.h"
+#include "py/runtime.h"
+#include "py/misc.h"
+#include "py/objtuple.h"
+
+#include "../../ndarray.h"
+#include "../../ulab.h"
+#include "../../ulab_tools.h"
+#include "integrate.h"
+
+ULAB_DEFINE_FLOAT_CONST(etolerance, MICROPY_FLOAT_CONST(1e-14), 0x283424dcUL, 0x3e901b2b29a4692bULL);
+
+static mp_float_t integrate_python_call(const mp_obj_type_t *type, mp_obj_t fun, mp_float_t x, mp_obj_t *fargs, uint8_t nparams) {
+    // Helper function for calculating the value of f(x, a, b, c, ...),
+    // where f is defined in python. Takes a float, returns a float.
+    // The array of mp_obj_t type must be supplied, as must the number of parameters (a, b, c...) in nparams
+    fargs[0] = mp_obj_new_float(x);
+    return mp_obj_get_float(MP_OBJ_TYPE_GET_SLOT(type, call)(fun, nparams+1, 0, fargs));
+}
+
+// sign helper function
+int sign(mp_float_t x) {
+    if (x >= 0) 
+        return 1;
+    else
+        return -1;
+}        
+
+// Tanh-Sinh, Sinh-Sinh and Exp-Sinh quadrature
+// https://www.genivia.com/qthsh.html
+
+// return optimized Exp-Sinh integral split point d
+mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t eps, mp_float_t d) {
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    mp_obj_t fargs[1];
+    mp_float_t h2 = integrate_python_call(type, fun, a + d/2, fargs, 0) - integrate_python_call(type, fun, (a + d*2)*4, fargs, 0);
+    int i = 1, j = 32;                   // j=32 is optimal to find r
+    if (MICROPY_FLOAT_C_FUN(isfinite)(h2) && MICROPY_FLOAT_C_FUN(fabs)(h2) > 1e-5) {    // if |h2| > 2^-16
+        mp_float_t r, fl, fr, h, s = 0, lfl, lfr, lr = 2;
+        do {                                  // find max j such that fl and fr are finite
+            j /= 2;
+            r = 1 << (i + j);
+            fl = integrate_python_call(type, fun, a + d/r, fargs, 0);
+            fr = integrate_python_call(type, fun, (a + d*r)*r*r, fargs, 0);
+            h = fl - fr;
+        } while (j > 1 && !MICROPY_FLOAT_C_FUN(isfinite)(h));
+        if (j > 1 && MICROPY_FLOAT_C_FUN(isfinite)(h) && sign(h) != sign(h2)) {
+            lfl = fl;                         // last fl=f(a+d/r)
+            lfr = fr;                         // last fr=f(a+d*r)*r*r
+            do {                              // bisect in 4 iterations
+                j /= 2;
+                r = 1 << (i + j);
+                fl = integrate_python_call(type, fun, a + d/r, fargs, 0);
+                fr = integrate_python_call(type, fun, (a + d*r)*r*r, fargs, 0);
+                h = fl - fr;
+                if (MICROPY_FLOAT_C_FUN(isfinite)(h)) {
+                    s += MICROPY_FLOAT_C_FUN(fabs)(h);  // sum |h| to remove noisy cases
+                    if (sign(h) == sign(h2)) {
+                        i += j;               // search right half
+                    }
+                    else {                    // search left half
+                        lfl = fl;             // record last fl=f(a+d/r)
+                        lfr = fr;             // record last fl=f(a+d*r)*r*r
+                        lr = r;               // record last r
+                    }
+                }
+            } while (j > 1);
+            if (s > eps) {                    // if sum of |h| > eps
+                h = lfl - lfr;                // use last fl and fr before the sign change
+                r = lr;                       // use last r before the sign change
+                if (h != 0)                   // if last diff != 0, back up r by one step
+                    r /= 2;
+                if (MICROPY_FLOAT_C_FUN(fabs)(lfl) < MICROPY_FLOAT_C_FUN(fabs)(lfr))
+                    d /= r;                   // move d closer to the finite endpoint
+                else
+                    d *= r;                   // move d closer to the infinite endpoint
+            }
+        }
+    }
+    return d;
+}
+
+
+// integrate function f, range a..b, max levels n, error tolerance eps
+mp_float_t quad(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint16_t n, mp_float_t eps, mp_float_t *e) {
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    mp_obj_t fargs[1];
+    const mp_float_t tol = 10*eps;
+    mp_float_t c = 0, d = 1, s, sign = 1, v, h = 2;
+    int k = 0, mode = 0;                   // Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
+    if (b < a) {                                // swap bounds
+        v = b;
+        b = a;
+        a = v;
+        sign = -1;
+    }
+    if (MICROPY_FLOAT_C_FUN(isfinite)(a) && MICROPY_FLOAT_C_FUN(isfinite)(b)) {
+        c = (a+b)/2;
+        d = (b-a)/2;
+        v = c;
+    }
+    else if (MICROPY_FLOAT_C_FUN(isfinite)(a)) {
+        mode = 1;                               // Exp-Sinh
+        d = exp_sinh_opt_d(fun, a, eps, d);
+        c = a;
+        v = a+d;
+    }
+    else if (MICROPY_FLOAT_C_FUN(isfinite)(b)) {
+        mode = 1;                               // Exp-Sinh
+        // d = -d;
+        d = exp_sinh_opt_d(fun, b, eps, -d);
+        sign = -sign;
+        c = b;
+        v = b+d;
+    }
+    else {
+        mode = 2;                               // Sinh-Sinh
+        v = 0;
+    }
+    s = integrate_python_call(type, fun, v, fargs, 0);
+    do {
+        mp_float_t p = 0, q, fp = 0, fm = 0, t, eh;
+        h /= 2;
+        t = eh = MICROPY_FLOAT_C_FUN(exp)(h);
+        if (k > 0)
+           eh *= eh;
+        if (mode == 0) {                        // Tanh-Sinh
+            do {
+                mp_float_t u = MICROPY_FLOAT_C_FUN(exp)(1/t-t);      // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
+                mp_float_t r = 2*u/(1+u);       // = 1 - tanh(sinh(j*h))
+                mp_float_t w = (t+1/t)*r/(1+u); // = cosh(j*h)/cosh(sinh(j*h))^2
+                mp_float_t x = d*r;
+                if (a+x > a) {                  // if too close to a then reuse previous fp
+                    mp_float_t y = integrate_python_call(type, fun, a+x, fargs, 0);
+                    if (MICROPY_FLOAT_C_FUN(isfinite)(y))
+                        fp = y;                 // if f(x) is finite, add to local sum
+                }
+                if (b-x < b) {                  // if too close to a then reuse previous fp
+                    mp_float_t y = integrate_python_call(type, fun, b-x, fargs, 0);
+                    if (MICROPY_FLOAT_C_FUN(isfinite)(y))
+                        fm = y;                 // if f(x) is finite, add to local sum
+                }
+                q = w*(fp+fm);
+                p += q;
+                t *= eh;
+            } while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p));
+        }
+        else {
+            t /= 2;
+            do {
+                mp_float_t r = MICROPY_FLOAT_C_FUN(exp)(t-.25/t);    // = exp(sinh(j*h))
+                mp_float_t x, y, w = r;
+                q = 0;
+                if (mode == 1) {                // Exp-Sinh
+                    x = c + d/r;
+                    if (x == c)                 // if x hit the finite endpoint then break
+                        break;
+                    y = integrate_python_call(type, fun, x, fargs, 0);
+                    if (MICROPY_FLOAT_C_FUN(isfinite)(y))    // if f(x) is finite, add to local sum
+                        q += y/w;
+                }
+                else {                          // Sinh-Sinh
+                    r = (r-1/r)/2;              // = sinh(sinh(j*h))
+                    w = (w+1/w)/2;              // = cosh(sinh(j*h))
+                    x = c - d*r;
+                    y = integrate_python_call(type, fun, x, fargs, 0);
+                    if (MICROPY_FLOAT_C_FUN(isfinite)(y))    // if f(x) is finite, add to local sum
+                        q += y*w;
+                }
+                x = c + d*r;
+                y = integrate_python_call(type, fun, x, fargs, 0);
+                if (MICROPY_FLOAT_C_FUN(isfinite)(y))        // if f(x) is finite, add to local sum
+                    q += y*w;
+                q *= t+.25/t;                   // q *= cosh(j*h)
+                p += q;
+                t *= eh;
+            } while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p));
+        }
+        v = s-p;
+        s += p;
+        ++k;
+    } while (MICROPY_FLOAT_C_FUN(fabs)(v) > tol*MICROPY_FLOAT_C_FUN(fabs)(s) && k <= n);
+    // return the error estimate by reference
+    *e = MICROPY_FLOAT_C_FUN(fabs)(v)/(MICROPY_FLOAT_C_FUN(fabs)(s)+eps);
+    return sign*d*s*h;                          // result with estimated relative error e
+}
+
+//| def quad(
+//|     fun: Callable[[float], float],
+//|     a: float,
+//|     b: float,
+//|     *,
+//|     levels: int = 6
+//|     eps: float = 1e-14,
+//| ) -> float:
+//|     """
+//|     :param callable f: The function to integrate
+//|     :param float a: The left side of the interval
+//|     :param float b: The right side of the interval
+//|     :param float levels: The number of levels to perform (6..7 is optimal)
+//|     :param float eps: The error tolerance value
+//|
+//|     Find a solution (zero) of the function ``f(x)`` on the interval
+//|     (``a``..``b``) using the bisection method.  The result is accurate to within
+//|     ``xtol`` unless more than ``maxiter`` steps are required."""
+//|     ...
+//|
+
+
+static mp_obj_t integrate_quad(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_levels, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 6} },
+        { MP_QSTR_eps, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = ULAB_REFERENCE_FLOAT_CONST(etolerance)} },
+    };
+
+    mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+
+    mp_obj_t fun = args[0].u_obj;
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    if(!MP_OBJ_TYPE_HAS_SLOT(type, call)) {
+        mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a function"));
+    }
+
+    mp_float_t a = mp_obj_get_float(args[1].u_obj);
+    mp_float_t b = mp_obj_get_float(args[2].u_obj);
+    uint16_t n = (uint16_t)args[3].u_int;
+    if(n < 0) {
+        mp_raise_ValueError(MP_ERROR_TEXT("levels should be > 0"));
+    }
+    mp_float_t eps = mp_obj_get_float(args[4].u_obj);
+    
+    mp_obj_t res[2];
+    mp_float_t e;
+    res[0] = mp_obj_new_float(quad(fun, a, b, n, eps, &e));
+    res[1] = mp_obj_new_float(e);
+    return mp_obj_new_tuple(2, res); 
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(integrate_quad_obj, 2, integrate_quad);
+
+
+// Romberg quadrature
+// This function is deprecated as of SciPy 1.12.0 and will be removed in SciPy 1.15.0. Please use scipy.integrate.quad instead. 
+// https://en.wikipedia.org/wiki/Romberg%27s_method, https://www.genivia.com/qthsh.html, 
+// https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romberg.html (which is different 
+// insofar it expects an array of function values). 
+
+mp_float_t qromb(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint16_t n, mp_float_t eps) {
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    mp_obj_t fargs[1];
+    mp_float_t R1[n], R2[n];
+    mp_float_t *Ro = &R1[0], *Ru = &R2[0];
+    mp_float_t h = b-a;
+    uint16_t i, j;
+    Ro[0] = (integrate_python_call(type, fun, a, fargs, 0) + integrate_python_call(type, fun, b, fargs, 0)) * h/2;
+    for (i = 1; i < n; ++i) {
+        unsigned long long k = 1UL << i;
+        unsigned long long s = 1;
+        mp_float_t sum = 0;
+        mp_float_t *Rt;
+        h /= 2;
+        for (j = 1; j < k; j += 2)
+            sum += integrate_python_call(type, fun, a+j*h, fargs, 0);
+        Ru[0] = h*sum + Ro[0]/2;
+        for (j = 1; j <= i; ++j) {
+            s <<= 2;
+            Ru[j] = (s*Ru[j-1] - Ro[j-1])/(s-1);
+        }
+        if (i > 2 && MICROPY_FLOAT_C_FUN(fabs)(Ro[i-1]-Ru[i]) <= eps*MICROPY_FLOAT_C_FUN(fabs)(Ru[i])+eps)
+            return Ru[i];
+        Rt = Ro;
+        Ro = Ru;
+        Ru = Rt;
+    }
+    return Ro[n-1];
+}
+
+//| def romberg(
+//|     fun: Callable[[float], float],
+//|     a: float,
+//|     b: float,
+//|     *,
+//|     steps: int = 100
+//|     eps: float = 1e-14,
+//| ) -> float:
+//|     """
+//|     :param callable f: The function to integrate
+//|     :param float a: The left side of the interval
+//|     :param float b: The right side of the interval
+//|     :param float steps: The number of equidistant steps
+//|     :param float eps: The tolerance value
+//|
+//|     Find a quadrature of the function ``f(x)`` on the interval
+//|     (``a``..``b``) using the Romberg method.  The result is accurate to within
+//|     ``eps`` unless more than ``steps`` steps are required."""
+//|     ...
+//|
+
+static mp_obj_t integrate_romberg(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_steps, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 100} },
+        { MP_QSTR_eps, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = ULAB_REFERENCE_FLOAT_CONST(etolerance)} },
+    };
+
+    mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+
+    mp_obj_t fun = args[0].u_obj;
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    if(!MP_OBJ_TYPE_HAS_SLOT(type, call)) {
+        mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a function"));
+    }
+
+    mp_float_t a = mp_obj_get_float(args[1].u_obj);
+    mp_float_t b = mp_obj_get_float(args[2].u_obj);
+    uint16_t steps = (uint16_t)args[3].u_int;
+    if(steps < 0) {
+        mp_raise_ValueError(MP_ERROR_TEXT("steps should be > 0"));
+    }
+    mp_float_t eps = mp_obj_get_float(args[4].u_obj);
+    
+    return mp_obj_new_float(qromb(fun, a, b, steps, eps)); 
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(integrate_romberg_obj, 2, integrate_romberg);
+
+
+// Adaptive Simpson quadrature
+// https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method, https://www.genivia.com/qthsh.html
+
+mp_float_t as(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, mp_float_t fa, mp_float_t fm,
+              mp_float_t fb, mp_float_t v, mp_float_t eps, int n, mp_float_t t) {
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    mp_obj_t fargs[1];
+    mp_float_t h = (b-a)/2;
+    mp_float_t f1 = integrate_python_call(type, fun, a + h/2, fargs, 0);
+    mp_float_t f2 = integrate_python_call(type, fun, b - h/2, fargs, 0);
+    mp_float_t sl = h*(fa + 4*f1 + fm)/6;
+    mp_float_t sr = h*(fm + 4*f2 + fb)/6;
+    mp_float_t s = sl+sr;
+    mp_float_t d = (s-v)/15;
+    mp_float_t m = a+h;
+    if (n <= 0 || MICROPY_FLOAT_C_FUN(fabs)(d) < eps)
+        return t + s + d; // note: fabs(d) can be used as error estimate
+    eps /= 2;
+    --n;
+    t = as(fun, a, m, fa, f1, fm, sl, eps, n, t);
+    return as(fun, m, b, fm, f2, fb, sr, eps, n, t);
+}
+
+mp_float_t qasi(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n, mp_float_t eps) {
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    mp_obj_t fargs[1];
+    mp_float_t fa = integrate_python_call(type, fun, a, fargs, 0);
+    mp_float_t fm = integrate_python_call(type, fun, (a+b)/2, fargs, 0);
+    mp_float_t fb = integrate_python_call(type, fun, b, fargs, 0);
+    mp_float_t v = (fa+4*fm+fb)*(b-a)/6;
+    return as(fun, a, b, fa, fm, fb, v, eps, n, 0);
+}
+
+//| def simpson(
+//|     fun: Callable[[float], float],
+//|     a: float,
+//|     b: float,
+//|     *,
+//|     steps: int = 100
+//|     eps: float = 1e-14,
+//| ) -> float:
+//|     """
+//|     :param callable f: The function to integrate
+//|     :param float a: The left side of the interval
+//|     :param float b: The right side of the interval
+//|     :param float steps: The number of equidistant steps
+//|     :param float eps: The tolerance value
+//|
+//|     Find a quadrature of the function ``f(x)`` on the interval
+//|     (``a``..``b``) using the Adaptive Simpson's method.  The result is accurate to within
+//|     ``eps`` unless more than ``steps`` steps are required."""
+//|     ...
+//|
+
+static mp_obj_t integrate_simpson(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_steps, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 100} },
+        { MP_QSTR_eps, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = ULAB_REFERENCE_FLOAT_CONST(etolerance)} },
+    };
+
+    mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+
+    mp_obj_t fun = args[0].u_obj;
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    if(!MP_OBJ_TYPE_HAS_SLOT(type, call)) {
+        mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a function"));
+    }
+
+    mp_float_t a = mp_obj_get_float(args[1].u_obj);
+    mp_float_t b = mp_obj_get_float(args[2].u_obj);
+    uint16_t steps = (uint16_t)args[3].u_int;
+    if(steps < 0) {
+        mp_raise_ValueError(MP_ERROR_TEXT("steps should be > 0"));
+    }
+    mp_float_t eps = mp_obj_get_float(args[4].u_obj);
+    
+    return mp_obj_new_float(qasi(fun, a, b, steps, eps)); 
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(integrate_simpson_obj, 2, integrate_simpson);
+
+
+// Adaptive Gauss-Kronrod (G10,K21) quadrature
+// https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula, https://www.genivia.com/qthsh.html
+
+mp_float_t gk(mp_float_t (*fun)(mp_float_t), mp_float_t c, mp_float_t d, mp_float_t *err) {
+// abscissas and weights pre-calculated with Legendre Stieltjes polynomials
+    static const mp_float_t abscissas[21] = {
+        0.00000000000000000e+00,
+        7.65265211334973338e-02,
+        1.52605465240922676e-01,
+        2.27785851141645078e-01,
+        3.01627868114913004e-01,
+        3.73706088715419561e-01,
+        4.43593175238725103e-01,
+        5.10867001950827098e-01,
+        5.75140446819710315e-01,
+        6.36053680726515025e-01,
+        6.93237656334751385e-01,
+        7.46331906460150793e-01,
+        7.95041428837551198e-01,
+        8.39116971822218823e-01,
+        8.78276811252281976e-01,
+        9.12234428251325906e-01,
+        9.40822633831754754e-01,
+        9.63971927277913791e-01,
+        9.81507877450250259e-01,
+        9.93128599185094925e-01,
+        9.98859031588277664e-01,
+    };
+    static const mp_float_t weights[21] = {
+        7.66007119179996564e-02,
+        7.63778676720807367e-02,
+        7.57044976845566747e-02,
+        7.45828754004991890e-02,
+        7.30306903327866675e-02,
+        7.10544235534440683e-02,
+        6.86486729285216193e-02,
+        6.58345971336184221e-02,
+        6.26532375547811680e-02,
+        5.91114008806395724e-02,
+        5.51951053482859947e-02,
+        5.09445739237286919e-02,
+        4.64348218674976747e-02,
+        4.16688733279736863e-02,
+        3.66001697582007980e-02,
+        3.12873067770327990e-02,
+        2.58821336049511588e-02,
+        2.03883734612665236e-02,
+        1.46261692569712530e-02,
+        8.60026985564294220e-03,
+        3.07358371852053150e-03,
+    };
+    static const mp_float_t gauss_weights[10] = {
+        1.52753387130725851e-01,
+        1.49172986472603747e-01,
+        1.42096109318382051e-01,
+        1.31688638449176627e-01,
+        1.18194531961518417e-01,
+        1.01930119817240435e-01,
+        8.32767415767047487e-02,
+        6.26720483341090636e-02,
+        4.06014298003869413e-02,
+        1.76140071391521183e-02,
+    };
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    mp_obj_t fargs[1];
+    mp_float_t p = 0; // kronrod quadrature sum
+    mp_float_t q = 0; // gauss quadrature sum
+    mp_float_t fp, fm;
+    mp_float_t e;
+    int i;
+    fp = integrate_python_call(type, fun, c, fargs, 0);
+    p = fp * weights[0];
+    for (i = 1; i < 21; i += 2) {
+        fp = integrate_python_call(type, fun, c + d * abscissas[i], fargs, 0);
+        fm = integrate_python_call(type, fun, c - d * abscissas[i], fargs, 0);
+        p += (fp + fm) * weights[i];
+        q += (fp + fm) * gauss_weights[i/2];
+    }
+    for (i = 2; i < 21; i += 2) {
+        fp = integrate_python_call(type, fun, c + d * abscissas[i], fargs, 0);
+        fm = integrate_python_call(type, fun, c - d * abscissas[i], fargs, 0);
+        p += (fp + fm) * weights[i];
+    }
+    *err = MICROPY_FLOAT_C_FUN(fabs)(p - q);
+    e = MICROPY_FLOAT_C_FUN(fabs)(2*p*1e-17); // optional, to take 1e-17 MachEps prec. into account
+    if (*err < e)
+        *err = e;
+    return p;
+}
+
+mp_float_t qakro(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n, mp_float_t tol, mp_float_t eps, mp_float_t *err) {
+    mp_float_t c = (a+b)/2;
+    mp_float_t d = (b-a)/2;
+    mp_float_t e;
+    mp_float_t r = gk(fun, c, d, &e);
+    mp_float_t s = d*r;
+    mp_float_t t = MICROPY_FLOAT_C_FUN(fabs)(s*tol);
+    if (tol == 0)
+        tol = t;
+    if (n > 0 && t < e && tol < e) {
+        s = qakro(fun, a, c, n-1, t/2, eps, err);
+        s += qakro(fun, c, b, n-1, t/2, eps, &e);
+        *err += e;
+        return s;
+    }
+    *err = e;
+    return s;
+}
+
+
+//| def quadgk(
+//|     fun: Callable[[float], float],
+//|     a: float,
+//|     b: float,
+//|     *,
+//|     order: int = 5
+//|     eps: float = 1e-14,
+//| ) -> float:
+//|     """
+//|     :param callable f: The function to integrate
+//|     :param float a: The left side of the interval
+//|     :param float b: The right side of the interval
+//|     :param float order: Order of quadrature integration. Default is 5.
+//|     :param float eps: The tolerance value
+//|
+//|     Find a quadrature of the function ``f(x)`` on the interval
+//|     (``a``..``b``) using the Adaptive Gauss-Kronrod method.  The result is accurate to within
+//|     ``eps`` unless a higher order than ``order`` is required."""
+//|     ...
+//|
+
+static mp_obj_t integrate_quadgk(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } },
+        { MP_QSTR_order, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 5} },
+        { MP_QSTR_eps, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = ULAB_REFERENCE_FLOAT_CONST(etolerance)} },
+    };
+
+    mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+
+    mp_obj_t fun = args[0].u_obj;
+    const mp_obj_type_t *type = mp_obj_get_type(fun);
+    if(!MP_OBJ_TYPE_HAS_SLOT(type, call)) {
+        mp_raise_TypeError(MP_ERROR_TEXT("first argument must be a function"));
+    }
+
+    mp_float_t a = mp_obj_get_float(args[1].u_obj);
+    mp_float_t b = mp_obj_get_float(args[2].u_obj);
+    uint16_t order = (uint16_t)args[3].u_int;
+    if(order < 0) {
+        mp_raise_ValueError(MP_ERROR_TEXT("levels should be > 0"));
+    }
+    mp_float_t eps = mp_obj_get_float(args[4].u_obj);
+    
+    mp_obj_t res[2];
+    mp_float_t e;
+    res[0] = mp_obj_new_float(qakro(fun, a, b, order, 0, eps, &e));
+    res[1] = mp_obj_new_float(e);
+    return mp_obj_new_tuple(2, res); 
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(integrate_quadgk_obj, 2, integrate_quadgk);
+
+static const mp_rom_map_elem_t ulab_scipy_integrate_globals_table[] = {
+    { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_integrate) },
+    { MP_ROM_QSTR(MP_QSTR_quad), MP_ROM_PTR(&integrate_quad_obj) },
+    { MP_ROM_QSTR(MP_QSTR_romberg), MP_ROM_PTR(&integrate_romberg_obj) },
+    { MP_ROM_QSTR(MP_QSTR_simpson), MP_ROM_PTR(&integrate_simpson_obj) },
+    { MP_ROM_QSTR(MP_QSTR_quadgk), MP_ROM_PTR(&integrate_quadgk_obj) },
+};
+
+static MP_DEFINE_CONST_DICT(mp_module_ulab_scipy_integrate_globals, ulab_scipy_integrate_globals_table);
+
+const mp_obj_module_t ulab_scipy_integrate_module = {
+    .base = { &mp_type_module },
+    .globals = (mp_obj_dict_t*)&mp_module_ulab_scipy_integrate_globals,
+};
+#if CIRCUITPY_ULAB
+MP_REGISTER_MODULE(MP_QSTR_ulab_dot_scipy_dot_integrate, ulab_scipy_integrate_module);
+#endif
+
diff -urN lvgl_micropython.O/ext_mod/ulab/code/scipy/integrate/integrate.h lvgl_micropython/ext_mod/ulab/code/scipy/integrate/integrate.h
--- lvgl_micropython.O/ext_mod/ulab/code/scipy/integrate/integrate.h	1970-01-01 01:00:00.000000000 +0100
+++ lvgl_micropython/ext_mod/ulab/code/scipy/integrate/integrate.h	2024-03-27 16:06:08.930794118 +0100
@@ -0,0 +1,44 @@
+
+/*
+ * This file is not part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2024 Harald Milz <hm@seneca.muc.de>
+ *
+*/
+
+#ifndef _SCIPY_INTEGRATE_
+#define _SCIPY_INTEGRATE_
+
+#include "../../ulab_tools.h"
+
+/*
+#ifndef     OPTIMIZE_EPSILON
+#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
+#define     OPTIMIZE_EPSILON      MICROPY_FLOAT_CONST(1.2e-7)
+#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
+#define     OPTIMIZE_EPSILON      MICROPY_FLOAT_CONST(2.3e-16)
+#endif
+#endif
+
+#define     OPTIMIZE_EPS          MICROPY_FLOAT_CONST(1.0e-4)
+#define     OPTIMIZE_NONZDELTA    MICROPY_FLOAT_CONST(0.05)
+#define     OPTIMIZE_ZDELTA       MICROPY_FLOAT_CONST(0.00025)
+#define     OPTIMIZE_ALPHA        MICROPY_FLOAT_CONST(1.0)
+#define     OPTIMIZE_BETA         MICROPY_FLOAT_CONST(2.0)
+#define     OPTIMIZE_GAMMA        MICROPY_FLOAT_CONST(0.5)
+#define     OPTIMIZE_DELTA        MICROPY_FLOAT_CONST(0.5)
+*/
+
+extern const mp_obj_module_t ulab_scipy_integrate_module;
+
+MP_DECLARE_CONST_FUN_OBJ_KW(optimize_quad_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(optimize_romberg_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(optimize_simpson_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(optimize_quadgk_obj);
+
+#endif /* _SCIPY_INTEGRATE_ */
+
diff -urN lvgl_micropython.O/ext_mod/ulab/code/scipy/scipy.c lvgl_micropython/ext_mod/ulab/code/scipy/scipy.c
--- lvgl_micropython.O/ext_mod/ulab/code/scipy/scipy.c	2024-11-25 14:59:24.500963171 +0100
+++ lvgl_micropython/ext_mod/ulab/code/scipy/scipy.c	2024-11-25 14:23:56.333974495 +0100
@@ -20,6 +20,8 @@
 #include "signal/signal.h"
 #include "special/special.h"
 #include "linalg/linalg.h"
+#include "integrate/integrate.h"
+
 
 #if ULAB_HAS_SCIPY
 
@@ -28,6 +30,9 @@
 
 static const mp_rom_map_elem_t ulab_scipy_globals_table[] = {
     { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_scipy) },
+    #if ULAB_SCIPY_HAS_INTEGRATE_MODULE
+        { MP_ROM_QSTR(MP_QSTR_integrate), MP_ROM_PTR(&ulab_scipy_integrate_module) },
+    #endif
     #if ULAB_SCIPY_HAS_LINALG_MODULE
         { MP_ROM_QSTR(MP_QSTR_linalg), MP_ROM_PTR(&ulab_scipy_linalg_module) },
     #endif
diff -urN lvgl_micropython.O/ext_mod/ulab/code/ulab.h lvgl_micropython/ext_mod/ulab/code/ulab.h
--- lvgl_micropython.O/ext_mod/ulab/code/ulab.h	2024-11-25 14:59:24.500963171 +0100
+++ lvgl_micropython/ext_mod/ulab/code/ulab.h	2024-11-25 14:25:05.789964308 +0100
@@ -398,6 +398,11 @@
 #define ULAB_NUMPY_HAS_WHERE            (1)
 #endif
 
+// the integrate module
+#ifndef ULAB_SCIPY_HAS_INTEGRATE_MODULE
+#define ULAB_SCIPY_HAS_INTEGRATE_MODULE        (1)
+#endif
+
 // the linalg module; functions of the linalg module still have
 // to be defined separately
 #ifndef ULAB_NUMPY_HAS_LINALG_MODULE
--- lvgl_micropython.O/lib/micropython/ports/esp32/mpconfigport.h	2024-11-25 15:05:25.362231890 +0100
+++ lvgl_micropython/lib/micropython/ports/esp32/mpconfigport.h	2024-11-25 14:50:30.967082119 +0100
@@ -65,7 +65,7 @@
 #define MICROPY_LONGINT_IMPL                (MICROPY_LONGINT_IMPL_MPZ)
 #define MICROPY_ERROR_REPORTING             (MICROPY_ERROR_REPORTING_NORMAL)
 #define MICROPY_WARNINGS                    (1)
-#define MICROPY_FLOAT_IMPL                  (MICROPY_FLOAT_IMPL_FLOAT)
+#define MICROPY_FLOAT_IMPL                  (MICROPY_FLOAT_IMPL_DOUBLE)
 #define MICROPY_STREAMS_POSIX_API           (1)
 #define MICROPY_USE_INTERNAL_ERRNO          (0) // errno.h from xtensa-esp32-elf/sys-include/sys
 #define MICROPY_USE_INTERNAL_PRINTF         (0) // ESP32 SDK requires its own printf
diff -urN lvgl_micropython.O/lib/micropython/ports/esp32/boards/ESP32_GENERIC_S3/mpconfigboard.h lvgl_micropython/lib/micropython/ports/esp32/boards/ESP32_GENERIC_S3/mpconfigboard.h
--- lvgl_micropython.O/lib/micropython/ports/esp32/boards/ESP32_GENERIC_S3/mpconfigboard.h	2024-11-25 14:57:36.520643625 +0100
+++ lvgl_micropython/lib/micropython/ports/esp32/boards/ESP32_GENERIC_S3/mpconfigboard.h	2024-11-25 14:42:31.576950765 +0100
@@ -9,3 +9,65 @@
 
 #define MICROPY_HW_I2C0_SCL                 (9)
 #define MICROPY_HW_I2C0_SDA                 (8)
+
+
+/* local additions for double float and complex handling */
+
+/* this is already in lib/micropython/ports/esp32/mpconfigport.h
+#ifdef MICROPY_CONFIG_ROM_LEVEL
+#undef MICROPY_CONFIG_ROM_LEVEL
+#endif
+#define MICROPY_CONFIG_ROM_LEVEL (MICROPY_CONFIG_ROM_LEVEL_EXTRA_FEATURES)
+
+#ifdef MICROPY_FLOAT_IMPL
+#undef MICROPY_FLOAT_IMPL
+#endif
+#define MICROPY_FLOAT_IMPL (MICROPY_FLOAT_IMPL_DOUBLE)
+*/
+
+#ifdef MICROPY_OBJ_REPR
+#undef MICROPY_OBJ_REPR
+#endif
+#define MICROPY_OBJ_REPR            (MICROPY_OBJ_REPR_A)
+
+// Whether to provide special math functions: math.{erf,erfc,gamma,lgamma}
+#ifdef MICROPY_PY_MATH_SPECIAL_FUNCTIONS
+#undef MICROPY_PY_MATH_SPECIAL_FUNCTIONS
+#endif
+#define MICROPY_PY_MATH_SPECIAL_FUNCTIONS (1)
+
+// Whether to provide math.factorial function
+#ifdef MICROPY_PY_MATH_FACTORIAL
+#undef MICROPY_PY_MATH_FACTORIAL
+#endif
+#define MICROPY_PY_MATH_FACTORIAL (1)
+
+// Whether to provide math.isclose function
+#ifdef MICROPY_PY_MATH_ISCLOSE
+#undef MICROPY_PY_MATH_ISCLOSE
+#endif
+#define MICROPY_PY_MATH_ISCLOSE (1)
+
+#ifdef MICROPY_PY_CMATH
+#undef MICROPY_PY_CMATH
+#endif
+#define MICROPY_PY_CMATH  (1)
+
+#ifdef ULAB_SUPPORTS_COMPLEX
+#undef ULAB_SUPPORTS_COMPLEX
+#endif
+#define ULAB_SUPPORTS_COMPLEX  (1)
+
+#ifdef MICROPY_OPT_MATH_FACTORIAL
+#undef MICROPY_OPT_MATH_FACTORIAL
+#endif
+#define MICROPY_OPT_MATH_FACTORIAL (1)
+
+#ifdef MICROPY_PY_BUILTINS_STR_UNICODE
+#undef MICROPY_PY_BUILTINS_STR_UNICODE
+#endif
+
+
+/* does not get resolved ...*/
+#define ULAB_NUMPY_HAS_LOG2             (1)
+