circuitpython-9.1.x.diff

diff -urN -x .git -x __pycache__ circuitpython/extmod/ulab/code/micropython.mk circuitpython.T/extmod/ulab/code/micropython.mk
--- circuitpython/extmod/ulab/code/micropython.mk	2024-03-24 17:06:55.801144545 +0100
+++ circuitpython.T/extmod/ulab/code/micropython.mk	2024-03-24 10:03:35.091506320 +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 -x .git -x __pycache__ circuitpython/extmod/ulab/code/scipy/integrate/integrate.c circuitpython.T/extmod/ulab/code/scipy/integrate/integrate.c
--- circuitpython/extmod/ulab/code/scipy/integrate/integrate.c	1970-01-01 01:00:00.000000000 +0100
+++ circuitpython.T/extmod/ulab/code/scipy/integrate/integrate.c	2024-03-24 16:54:54.253321641 +0100
@@ -0,0 +1,639 @@
+/*
+ * 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 -x .git -x __pycache__ circuitpython/extmod/ulab/code/scipy/integrate/integrate.h circuitpython.T/extmod/ulab/code/scipy/integrate/integrate.h
--- circuitpython/extmod/ulab/code/scipy/integrate/integrate.h	1970-01-01 01:00:00.000000000 +0100
+++ circuitpython.T/extmod/ulab/code/scipy/integrate/integrate.h	2024-03-24 15:45:45.517189057 +0100
@@ -0,0 +1,43 @@
+
+/*
+ * 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 -x .git -x __pycache__ circuitpython/extmod/ulab/code/scipy/scipy.c circuitpython.T/extmod/ulab/code/scipy/scipy.c
--- circuitpython/extmod/ulab/code/scipy/scipy.c	2024-03-24 17:06:55.805144419 +0100
+++ circuitpython.T/extmod/ulab/code/scipy/scipy.c	2024-03-24 10:11:48.370193665 +0100
@@ -20,6 +20,7 @@
 #include "signal/signal.h"
 #include "special/special.h"
 #include "linalg/linalg.h"
+#include "integrate/integrate.h"
 
 #if ULAB_HAS_SCIPY
 
@@ -28,6 +29,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 -x .git -x __pycache__ circuitpython/extmod/ulab/code/ulab.h circuitpython.T/extmod/ulab/code/ulab.h
--- circuitpython/extmod/ulab/code/ulab.h	2024-03-24 17:06:55.805144419 +0100
+++ circuitpython.T/extmod/ulab/code/ulab.h	2024-03-24 10:10:21.842160534 +0100
@@ -698,6 +698,10 @@
 #endif
 
 // scipy modules
+#ifndef ULAB_SCIPY_HAS_INTEGRATE_MODULE
+#define ULAB_SCIPY_HAS_INTEGRATE_MODULE        (1)
+#endif
+
 #ifndef ULAB_SCIPY_HAS_LINALG_MODULE
 #define ULAB_SCIPY_HAS_LINALG_MODULE        (1)
 #endif
diff -urN -x .git -x __pycache__ circuitpython/ports/atmel-samd/boards/feather_m4_express/mpconfigboard.h circuitpython.T/ports/atmel-samd/boards/feather_m4_express/mpconfigboard.h
--- circuitpython/ports/atmel-samd/boards/feather_m4_express/mpconfigboard.h	2024-03-24 17:06:05.706791128 +0100
+++ circuitpython.T/ports/atmel-samd/boards/feather_m4_express/mpconfigboard.h	2024-03-21 12:57:49.546596927 +0100
@@ -23,3 +23,67 @@
 // USB is always used internally so skip the pin objects for it.
 #define IGNORE_PIN_PA24     1
 #define IGNORE_PIN_PA25     1
+
+/* local additions for double float and complex handling */
+
+/*
+#ifndef MICROPY_CONFIG_ROM_LEVEL
+#define MICROPY_CONFIG_ROM_LEVEL (MICROPY_CONFIG_ROM_LEVEL_EXTRA_FEATURES)
+#endif
+*/
+
+#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             (0)
+
+
+#pragma GCC diagnostic ignored "-Woverflow"
+
+
diff -urN -x .git -x __pycache__ circuitpython/ports/atmel-samd/boards/feather_m4_express/mpconfigboard.mk circuitpython.T/ports/atmel-samd/boards/feather_m4_express/mpconfigboard.mk
--- circuitpython/ports/atmel-samd/boards/feather_m4_express/mpconfigboard.mk	2024-03-24 17:06:05.706791128 +0100
+++ circuitpython.T/ports/atmel-samd/boards/feather_m4_express/mpconfigboard.mk	2024-03-20 09:22:09.423869743 +0100
@@ -15,6 +15,25 @@
 CIRCUITPY_SYNTHIO = 0
 CIRCUITPY_JPEGIO = 0
 
+# some more to make room for math
+CIRCUITPY__BLEIO = 0
+CIRCUITPY_AUDIOBUSIO = 0
+CIRCUITPY_AUDIOCORE = 0
+CIRCUITPY_AUDIOIO = 0
+CIRCUITPY_AUDIOMIXER = 0
+CIRCUITPY_AUDIOMP3 = 0
+CIRCUITPY_GETPASS = 0
+CIRCUITPY_GIFIO = 0
+CIRCUITPY_KEYPAD = 0
+CIRCUITPY_ONEWIREIO = 0
+CIRCUITPY_PARALLELDISPLAY = 0
+CIRCUITPY_PS2IO = 0
+CIRCUITPY_ROTARYIO = 0
+CIRCUITPY_SHARPDISPLAY = 0
+CIRCUITPY_USB_MIDI = 0
+CIRCUITPY_ZLIB = 0
+
+
 # We don't have room for the fonts for terminalio for certain languages,
 # so turn off terminalio, and if it's off and displayio is on,
 # force a clean build.
diff -urN -x .git -x __pycache__ circuitpython/ports/atmel-samd/Makefile circuitpython.T/ports/atmel-samd/Makefile
--- circuitpython/ports/atmel-samd/Makefile	2024-03-24 17:06:05.698791404 +0100
+++ circuitpython.T/ports/atmel-samd/Makefile	2024-03-21 13:28:12.020574999 +0100
@@ -127,6 +127,10 @@
 
 CFLAGS += $(INC) -Wall -Werror -std=gnu11 -nostdlib -fshort-enums $(BASE_CFLAGS) $(CFLAGS_MOD) $(COPT) -Werror=missing-prototypes
 
+ifneq ($(INTERNAL_LIBM_DBL),1)
+CFLAGS += -Werror
+endif
+
 ifeq ($(CHIP_FAMILY), samd21)
 CFLAGS += \
 	-mthumb \
diff -urN -x .git -x __pycache__ circuitpython/ports/atmel-samd/mpconfigport.mk circuitpython.T/ports/atmel-samd/mpconfigport.mk
--- circuitpython/ports/atmel-samd/mpconfigport.mk	2024-03-24 17:06:06.066778725 +0100
+++ circuitpython.T/ports/atmel-samd/mpconfigport.mk	2024-03-20 09:22:09.423869743 +0100
@@ -2,7 +2,9 @@
 # parameters that vary based on chip and/or board.
 LD_TEMPLATE_FILE = boards/common.template.ld
 
-INTERNAL_LIBM = 1
+INTERNAL_LIBM=
+# INTERNAL_LIBM = 1
+# INTERNAL_LIBM_DBL = 1
 
 # Number of USB endpoint pairs.
 USB_NUM_ENDPOINT_PAIRS = 8
diff -urN -x .git -x __pycache__ circuitpython/ports/raspberrypi/boards/adafruit_feather_rp2040/mpconfigboard.h circuitpython.T/ports/raspberrypi/boards/adafruit_feather_rp2040/mpconfigboard.h
--- circuitpython/ports/raspberrypi/boards/adafruit_feather_rp2040/mpconfigboard.h	2024-03-24 17:06:05.834786717 +0100
+++ circuitpython.T/ports/raspberrypi/boards/adafruit_feather_rp2040/mpconfigboard.h	2024-03-21 13:49:50.777482266 +0100
@@ -12,3 +12,49 @@
 
 #define DEFAULT_UART_BUS_RX (&pin_GPIO1)
 #define DEFAULT_UART_BUS_TX (&pin_GPIO0)
+
+/* local additions for double float and complex handling */
+/*
+#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)
+
diff -urN -x .git -x __pycache__ circuitpython/ports/raspberrypi/Makefile circuitpython.T/ports/raspberrypi/Makefile
--- circuitpython/ports/raspberrypi/Makefile	2024-03-24 17:06:05.834786717 +0100
+++ circuitpython.T/ports/raspberrypi/Makefile	2024-03-21 13:27:32.796442507 +0100
@@ -144,6 +144,11 @@
         -I../../supervisor/shared/usb \
         -I$(BUILD)
 
+ifeq ($(INTERNAL_LIBM_DBL),1)
+INC += \
+        -isystem sdk/src/rp2_common/pico_double/include 
+endif
+
 # Pico specific configuration
 CFLAGS += -DRASPBERRYPI -DPICO_ON_DEVICE=1 -DPICO_NO_BINARY_INFO=0 -DPICO_TIME_DEFAULT_ALARM_POOL_DISABLED=0 -DPICO_DIVIDER_CALL_IDIV0=0 -DPICO_DIVIDER_CALL_LDIV0=0 -DPICO_DIVIDER_HARDWARE=1 -DPICO_DOUBLE_ROM=1 -DPICO_FLOAT_ROM=1 -DPICO_MULTICORE=1 -DPICO_BITS_IN_RAM=0 -DPICO_DIVIDER_IN_RAM=0 -DPICO_DOUBLE_PROPAGATE_NANS=0 -DPICO_DOUBLE_IN_RAM=0 -DPICO_MEM_IN_RAM=0 -DPICO_FLOAT_IN_RAM=0 -DPICO_FLOAT_PROPAGATE_NANS=1 -DPICO_NO_FLASH=0 -DPICO_COPY_TO_RAM=0 -DPICO_DISABLE_SHARED_IRQ_HANDLERS=0 -DPICO_NO_BI_BOOTSEL_VIA_DOUBLE_RESET=0 -DDVI_1BPP_BIT_REVERSE=0
 OPTIMIZATION_FLAGS ?= -O3
@@ -173,7 +178,11 @@
 # Remove -Wno-stringop-overflow after we can test with CI's GCC 10. Mac's looks weird.
 DISABLE_WARNINGS = -Wno-stringop-overflow -Wno-cast-align
 
-CFLAGS += $(INC) -Wall -Werror -std=gnu11 -fshort-enums $(BASE_CFLAGS) $(CFLAGS_MOD) $(COPT) $(DISABLE_WARNINGS) -Werror=missing-prototypes
+CFLAGS += $(INC) -Wall -std=gnu11 -fshort-enums $(BASE_CFLAGS) $(CFLAGS_MOD) $(COPT) $(DISABLE_WARNINGS) -Werror=missing-prototypes
+
+ifneq ($(INTERNAL_LIBM_DBL),1)
+CFLAGS += -Werror
+endif
 
 CFLAGS += \
 	-march=armv6-m \
@@ -225,14 +234,19 @@
 	src/rp2_common/pico_fix/rp2040_usb_device_enumeration/rp2040_usb_device_enumeration.c \
 	src/rp2_common/pico_float/float_init_rom.c \
 	src/rp2_common/pico_float/float_math.c \
-	src/rp2_common/pico_multicore/multicore.c \
+    src/rp2_common/pico_multicore/multicore.c \
 	src/rp2_common/pico_platform/platform.c \
 	src/rp2_common/pico_printf/printf.c \
 	src/rp2_common/pico_runtime/runtime.c \
 	src/rp2_common/pico_stdio/stdio.c \
 	src/rp2_common/pico_stdlib/stdlib.c \
 	src/rp2_common/pico_unique_id/unique_id.c \
-	$(SRC_SDK_CYW43) \
+	$(SRC_SDK_CYW43) \ 
+
+ifeq ($(INTERNAL_LIBM_DBL),1)
+SRC_SDK += \
+    src/rp2_common/pico_double/double_math.c
+endif
 
 SRC_SDK := $(addprefix sdk/, $(SRC_SDK))
 $(patsubst %.c,$(BUILD)/%.o,$(SRC_SDK) $(SRC_CYW43)): CFLAGS += -Wno-missing-prototypes -Wno-undef -Wno-unused-function -Wno-nested-externs -Wno-strict-prototypes -Wno-double-promotion -Wno-sign-compare -Wno-unused-variable -Wno-strict-overflow -Ilib/cyw43-driver
diff -urN -x .git -x __pycache__ circuitpython/ports/raspberrypi/mpconfigport.mk circuitpython.T/ports/raspberrypi/mpconfigport.mk
--- circuitpython/ports/raspberrypi/mpconfigport.mk	2024-03-24 17:06:05.854786028 +0100
+++ circuitpython.T/ports/raspberrypi/mpconfigport.mk	2024-03-20 09:22:09.423869743 +0100
@@ -46,7 +46,8 @@
 
 CIRCUITPY_AUDIOMIXER ?= 1
 
-INTERNAL_LIBM = 1
+INTERNAL_LIBM=
+#INTERNAL_LIBM_DBL = 1
 
 CIRCUITPY_BUILD_EXTENSIONS ?= uf2
 
diff -urN -x .git -x __pycache__ circuitpython/py/circuitpy_defns.mk circuitpython.T/py/circuitpy_defns.mk
--- circuitpython/py/circuitpy_defns.mk	2024-03-24 17:06:06.098777625 +0100
+++ circuitpython.T/py/circuitpy_defns.mk	2024-03-20 09:22:09.423869743 +0100
@@ -28,14 +28,11 @@
 # Common compile warnings.
 
 BASE_CFLAGS = \
-	-fsingle-precision-constant \
 	-fno-strict-aliasing \
-	-Wdouble-promotion \
 	-Wimplicit-fallthrough=2 \
 	-Wno-endif-labels \
 	-Wstrict-prototypes \
 	-Werror-implicit-function-declaration \
-	-Wfloat-equal \
 	-Wundef \
 	-Wshadow \
 	-Wwrite-strings \
@@ -58,6 +55,19 @@
 #        -ftime-report
 #        -H
 
+# only if we do not build _DOUBLE
+# and even then ... 
+# by the way this makes all these hardcoded f's redundant
+# they should all be thrown out the window. 
+ifdef INTERNAL_LIBM 
+ifndef INTERNAL_LIBM_DBL
+BASE_CFLAGS += \
+	-fsingle-precision-constant \
+	-Wdouble-promotion \
+	-Wfloat-equal 
+endif
+endif
+
 # Micropython's implementation of <string.h> routines is incompatible with
 # "fortify source", enabled by default on gentoo's crossdev arm-none-eabi-gcc
 # gcc version 12.3.1 20230526 (Gentoo 12.3.1_p20230526 p2). Unconditionally disable it.
@@ -874,6 +884,60 @@
 $(patsubst %.c,$(BUILD)/%.o,$(SRC_LIBM)): CFLAGS += -Wno-missing-prototypes
 endif
 
+# and if INTERNAL_LIBM_DBL == 1, we also include the libm/libm_dbl objects. 
+ifeq ($(INTERNAL_LIBM_DBL),1)
+SRC_LIBM += \
+$(addprefix lib/,\
+	libm_dbl/acos.c \
+	libm_dbl/acosh.c \
+	libm_dbl/asin.c \
+	libm_dbl/asinh.c \
+	libm_dbl/atan2.c \
+	libm_dbl/atan.c \
+	libm_dbl/atanh.c \
+	libm_dbl/ceil.c \
+	libm_dbl/copysign.c \
+	libm_dbl/__cos.c \
+	libm_dbl/cos.c \
+	libm_dbl/cosh.c \
+	libm_dbl/erf.c \
+	libm_dbl/exp.c \
+	libm_dbl/expm1.c \
+	libm_dbl/__expo2.c \
+	libm_dbl/floor.c \
+	libm_dbl/fmod.c \
+	libm_dbl/__fpclassify.c \
+	libm_dbl/frexp.c \
+	libm_dbl/ldexp.c \
+	libm_dbl/lgamma.c \
+	libm_dbl/log10.c \
+	libm_dbl/log1p.c \
+	libm_dbl/log.c \
+	libm_dbl/modf.c \
+	libm_dbl/nearbyint.c \
+	libm_dbl/pow.c \
+	libm_dbl/__rem_pio2.c \
+	libm_dbl/__rem_pio2_large.c \
+	libm_dbl/rint.c \
+	libm_dbl/round.c \
+	libm_dbl/scalbn.c \
+	libm_dbl/__signbit.c \
+	libm_dbl/__sin.c \
+	libm_dbl/sin.c \
+	libm_dbl/sinh.c \
+	libm_dbl/sqrt.c \
+	libm_dbl/__tan.c \
+	libm_dbl/tan.c \
+	libm_dbl/tanh.c \
+	libm_dbl/tgamma.c \
+	libm_dbl/trunc.c \
+	)
+# we include the special functions without ULAB. 
+endif # INTERNAL_LIBM == 2
+
+	# libm_dbl/thumb_vfp_sqrt.c \
+
+
 # Sources used in all ports except unix.
 SRC_CIRCUITPY_COMMON = \
 	shared/libc/string0.c \
diff -urN -x .git -x __pycache__ circuitpython/py/circuitpy_mpconfig.h circuitpython.T/py/circuitpy_mpconfig.h
--- circuitpython/py/circuitpy_mpconfig.h	2024-03-24 17:06:06.098777625 +0100
+++ circuitpython.T/py/circuitpy_mpconfig.h	2024-03-20 09:22:09.423869743 +0100
@@ -368,7 +368,7 @@
 #endif
 
 #ifndef ULAB_SUPPORTS_COMPLEX
-#define ULAB_SUPPORTS_COMPLEX (0)
+#define ULAB_SUPPORTS_COMPLEX (1)
 #endif
 
 #if CIRCUITPY_ULAB
diff -urN -x .git -x __pycache__ circuitpython/py/obj.h circuitpython.T/py/obj.h
--- circuitpython/py/obj.h	2024-03-24 17:06:05.886784925 +0100
+++ circuitpython.T/py/obj.h	2024-03-20 09:22:09.423869743 +0100
@@ -194,7 +194,9 @@
 
 static inline bool mp_obj_is_float(mp_const_obj_t o) {
     // Ensure that 32-bit arch can only use single precision.
+#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
     MP_STATIC_ASSERT(sizeof(mp_float_t) <= sizeof(mp_obj_t));
+#endif	
 
     return (((mp_uint_t)(o)) & 3) == 2 && (((mp_uint_t)(o)) & 0xff800007) != 0x00000006;
 }
diff -urN -x .git -x __pycache__ circuitpython/py/unicode.c circuitpython.T/py/unicode.c
--- circuitpython/py/unicode.c	2024-03-24 17:06:05.894784649 +0100
+++ circuitpython.T/py/unicode.c	2024-03-21 13:12:51.849485297 +0100
@@ -64,6 +64,18 @@
     AT_PR, AT_LX, AT_LX, AT_LX, AT_LX, AT_LX, AT_LX, AT_LO,
     AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO,
     AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO,
+    AT_LO, AT_LO, AT_LO, AT_PR, AT_PR, AT_PR, AT_PR, 0,
+	AT_SP, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR,
+    AT_PR, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR,
+    AT_DI, AT_DI, AT_DI, AT_DI, AT_DI, AT_DI, AT_DI, AT_DI,
+    AT_DI, AT_DI, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR,
+    AT_PR, AT_UX, AT_UX, AT_UX, AT_UX, AT_UX, AT_UX, AT_UP,
+    AT_UP, AT_UP, AT_UP, AT_UP, AT_UP, AT_UP, AT_UP, AT_UP,
+    AT_UP, AT_UP, AT_UP, AT_UP, AT_UP, AT_UP, AT_UP, AT_UP,
+    AT_UP, AT_UP, AT_UP, AT_PR, AT_PR, AT_PR, AT_PR, AT_PR,
+    AT_PR, AT_LX, AT_LX, AT_LX, AT_LX, AT_LX, AT_LX, AT_LO,
+    AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO,
+    AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO, AT_LO,
     AT_LO, AT_LO, AT_LO, AT_PR, AT_PR, AT_PR, AT_PR, 0
 };
 
@@ -117,41 +129,41 @@
 
 // Be aware: These unichar_is* functions are actually ASCII-only!
 bool unichar_isspace(unichar c) {
-    return c < 128 && (attr[c] & FL_SPACE) != 0;
+    return c < 224 && (attr[c] & FL_SPACE) != 0;
 }
 
 bool unichar_isalpha(unichar c) {
-    return c < 128 && (attr[c] & FL_ALPHA) != 0;
+    return c < 224 && (attr[c] & FL_ALPHA) != 0;
 }
 
 /* unused
 bool unichar_isprint(unichar c) {
-    return c < 128 && (attr[c] & FL_PRINT) != 0;
+    return c < 224 && (attr[c] & FL_PRINT) != 0;
 }
 */
 
 bool unichar_isdigit(unichar c) {
-    return c < 128 && (attr[c] & FL_DIGIT) != 0;
+    return c < 224 && (attr[c] & FL_DIGIT) != 0;
 }
 
 bool unichar_isxdigit(unichar c) {
-    return c < 128 && (attr[c] & FL_XDIGIT) != 0;
+    return c < 224 && (attr[c] & FL_XDIGIT) != 0;
 }
 
 bool unichar_isident(unichar c) {
-    return c < 128 && ((attr[c] & (FL_ALPHA | FL_DIGIT)) != 0 || c == '_');
+    return c < 224 && ((attr[c] & (FL_ALPHA | FL_DIGIT)) != 0 || c == '_');
 }
 
 bool unichar_isalnum(unichar c) {
-    return c < 128 && ((attr[c] & (FL_ALPHA | FL_DIGIT)) != 0);
+    return c < 224 && ((attr[c] & (FL_ALPHA | FL_DIGIT)) != 0);
 }
 
 bool unichar_isupper(unichar c) {
-    return c < 128 && (attr[c] & FL_UPPER) != 0;
+    return c < 224 && (attr[c] & FL_UPPER) != 0;
 }
 
 bool unichar_islower(unichar c) {
-    return c < 128 && (attr[c] & FL_LOWER) != 0;
+    return c < 224 && (attr[c] & FL_LOWER) != 0;
 }
 
 unichar unichar_tolower(unichar c) {
diff -urN -x .git -x __pycache__ circuitpython/shared-bindings/busdisplay/BusDisplay.c circuitpython.T/shared-bindings/busdisplay/BusDisplay.c
--- circuitpython/shared-bindings/busdisplay/BusDisplay.c	2024-03-24 17:06:06.098777625 +0100
+++ circuitpython.T/shared-bindings/busdisplay/BusDisplay.c	2024-03-20 09:20:33.123349901 +0100
@@ -104,7 +104,7 @@
 //|                            b"\x11\x80\x78"# Exit Sleep then delay 0x78 (120ms)
 //|                            b"\x29\x81\xaa\x78"# Display on then delay 0x78 (120ms)
 //|                           )
-//|            display = busdisplay.BusDisplay(display_bus, init_sequence, width=320, height=240)
+//|            display = displayio.Display(display_bus, init_sequence, width=320, height=240)
 //|
 //|         The first command is 0xe1 with 15 (0xf) parameters following. The second is 0x11 with 0
 //|         parameters and a 120ms (0x78) delay. The third command is 0x29 with one parameter 0xaa and a
diff -urN -x .git -x __pycache__ circuitpython/shared-bindings/rgbmatrix/RGBMatrix.c circuitpython.T/shared-bindings/rgbmatrix/RGBMatrix.c
--- circuitpython/shared-bindings/rgbmatrix/RGBMatrix.c	2024-03-24 17:06:05.914783960 +0100
+++ circuitpython.T/shared-bindings/rgbmatrix/RGBMatrix.c	2024-03-21 13:26:49.396295952 +0100
@@ -390,7 +390,11 @@
     common_hal_rgbmatrix_rgbmatrix_deinit(self_in);
 }
 
+#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
+static mp_float_t rgbmatrix_rgbmatrix_get_brightness_proto(mp_obj_t self_in) {
+#else
 static float rgbmatrix_rgbmatrix_get_brightness_proto(mp_obj_t self_in) {
+#endif
     return common_hal_rgbmatrix_rgbmatrix_get_paused(self_in) ? 0.0f : 1.0f;
 }
 
diff -urN -x .git -x __pycache__ circuitpython/shared-module/fontio/BuiltinFont.c circuitpython.T/shared-module/fontio/BuiltinFont.c
--- circuitpython/shared-module/fontio/BuiltinFont.c	2024-03-24 17:06:05.926783546 +0100
+++ circuitpython.T/shared-module/fontio/BuiltinFont.c	2024-03-21 13:18:51.274685867 +0100
@@ -43,7 +43,7 @@
 }
 
 uint8_t fontio_builtinfont_get_glyph_index(const fontio_builtinfont_t *self, mp_uint_t codepoint) {
-    if (codepoint >= 0x20 && codepoint <= 0x7e) {
+    if (codepoint >= 0x20 && codepoint <= 0xde) {
         return codepoint - 0x20;
     }
     // Do a linear search of the mapping for unicode.
@@ -53,7 +53,7 @@
         unichar potential_c = utf8_get_char(j);
         j = utf8_next_char(j);
         if (codepoint == potential_c) {
-            return 0x7f - 0x20 + k;
+            return 0xdf - 0x20 + k;
         }
         k++;
     }
diff -urN -x .git -x __pycache__ circuitpython/shared-module/terminalio/Terminal.c circuitpython.T/shared-module/terminalio/Terminal.c
--- circuitpython/shared-module/terminalio/Terminal.c	2024-03-24 17:06:05.930783408 +0100
+++ circuitpython.T/shared-module/terminalio/Terminal.c	2024-03-21 13:12:18.881375889 +0100
@@ -98,8 +98,8 @@
             continue;
         }
         // Always handle ASCII.
-        if (c < 128) {
-            if (c >= 0x20 && c <= 0x7e) {
+        if (c < 224) {
+            if (c >= 0x20 && c <= 0xde) {
                 uint8_t tile_index = fontio_builtinfont_get_glyph_index(self->font, c);
                 common_hal_displayio_tilegrid_set_tile(self->scroll_area, self->cursor_x, self->cursor_y, tile_index);
                 self->cursor_x++;
diff -urN -x .git -x __pycache__ circuitpython/tools/gen_display_resources.py circuitpython.T/tools/gen_display_resources.py
--- circuitpython/tools/gen_display_resources.py	2024-03-24 17:06:06.014780514 +0100
+++ circuitpython.T/tools/gen_display_resources.py	2024-03-20 09:22:09.423869743 +0100
@@ -286,7 +286,7 @@
     .in_group = true
 }};
 """.format(
-        len(all_characters), tile_x, tile_y
+        len(all_characters) * 2, tile_x, tile_y
     )
 )
 
@@ -316,7 +316,7 @@
     .in_group = true
 }};
 """.format(
-        len(all_characters), tile_x, tile_y
+        len(all_characters) * 2, tile_x, tile_y
     )
 )
 
@@ -324,14 +324,19 @@
     """\
 const uint32_t font_bitmap_data[{}] = {{
 """.format(
-        bytes_per_row * tile_y // 4
+        bytes_per_row * tile_y // 4 * 2
     )
 )
 
+# for each 18 words written, add another 18 words with inverted data. 
+extras = ""
 for i, word in enumerate(struct.iter_unpack(">I", b)):
     c_file.write("0x{:08x}, ".format(word[0]))
+    extras += "0x{:08x}, ".format(word[0] ^ 0xFFFFFFFF)
     if (i + 1) % (bytes_per_row // 4) == 0:
-        c_file.write("\n")
+        extras += "\n"
+        c_file.write(extras)
+        extras = ""
 
 c_file.write(
     """\
@@ -354,7 +359,7 @@
     .read_only = true
 }};
 """.format(
-        len(all_characters) * tile_x, tile_y, bytes_per_row / 4
+        2 * len(all_characters) * tile_x, tile_y, bytes_per_row / 4 * 2
     )
 )