github actions: Fix paths directives

I noticed on #122 that no actions were run.  I believe this is because
the paths directives had no wildcards, contrary to the examples at
https://help.github.com/en/actions/reference/workflow-syntax-for-github-actions#onpushpull_requestpaths

Eliminating the paths: restrictions entirely may be an equally valid
course of action.

.github/workflows/build.yml

@@ -3,6 +3,10 @@ name: Build CI
 on:
   push:
   pull_request:
+    paths:
+    - 'code/**'
+    - 'tests/**'
+    - '.github/workflows/**'
   release:
     types: [published]
   check_suite:

.gitignore

@@ -0,0 +1,3 @@
+/micropython
+/*.exp
+/*.out

README.md

@@ -1,25 +1,84 @@
-# micropython-ulab
+# ulab
 
-ulab is a numpy-like array manipulation library for micropython. 
-The module is written in C, defines compact containers for numerical 
-data, and is fast. 
+`ulab` is a `numpy`-like array manipulation library for `micropython` and `CircuitPython`.
+The module is written in C, defines compact containers for numerical
+data, and is fast. The library is a software-only standard `micropython` user module,
+i.e., it has no hardware dependencies, and can be compiled for any platform. 
+The `float` implementation of `micropython` (`float`, or `double`) is automatically detected.
 
-Documentation can be found under https://micropython-ulab.readthedocs.io/en/latest/
-The source for the manual is in https://github.com/v923z/micropython-ulab/blob/master/docs/ulab-manual.ipynb,
-while developer help is in https://github.com/v923z/micropython-ulab/blob/master/docs/ulab.ipynb.
+# Supported functions
+
+
+## ndarray
+
+`ulab` implements `numpy`'s `ndarray` with the `==`, `!=`, `<`, `<=`, `>`, `>=`, `+`, `-`, `/`, `*`, and `**` binary operators, and the `len`, `~`, `-`, `+`, `abs` unary operators that operate element-wise. Type-aware `ndarray`s can be initialised from any `micropython` iterable, lists of iterables, or by means of the `ones`, `zeros`, `eye`, or `linspace` functions. 
+
+`ndarray`s can be iterated on, and have a number of their own methods, such as `shape`, `reshape`, `transpose`, `size`, and `itemsize`.
+
+## Modules
+
+In addition to the `ndarray`'s operators and methods, seven modules define a great number of functions that can take `ndarray`s or `micropython` iterables as their arguments. If flash space is a concern, unnecessary sub-modules can be excluded from the compiled firmware with a pre-processor switch. 
+
+### vector
+
+The `vector` sub-module implements all functions of the `math` package (e.g., `acos`, `acosh`, ..., `tan`, `tanh`) of `micropython` for `ndarray`s and iterables. In addition, it also provided tools for vectorising generic, user-defined `python` functions. 
+
+### numerical
+
+The `numerical` sub-module defines the `roll`, `flip`, `diff`, `sort` and `argsort` functions for `ndarray`s, and, in addition, the `min`, `max`, `argmin`, `argmax`, `sum`, `mean`, `std` functions that work with `ndarray`s, as well as generic one-dimensional iterables.
+
+### linalg
+
+The `linalg` sub-module implements functions for matrix inversion, dot product, and the calculation of the determinant, eigenvalues, eigenvectors, Cholesky decomposition, and trace. 
+
+### poly
+
+The `poly` sub-module defines the `polyval`, and `polyfit` functions from `numpy`.
+
+### fft
+
+The `fft` sub-module implements the fast Fourier transform, and its inverse for one-dimensional `ndarray`s.
+
+### filter
+
+The `filter` sub-module implements one-dimensional convolution.
+
+### compare
+
+The `compare` sub-module contains the implementation of the `equal`, `not_equal`, `minimum`, `maximum`, and `clip` functions.
+
+### approx
+
+The `approx` sub-module contains the implementation of the `interp` function of `numpy`, and `newton`, `bisect`, and `fmin` from `scipy`.
+
+### extras
+
+The `extras` sub-module is meant as a user-extendable module, and currently implements the `spectrogram` function of `scipy`. 
+
+# Finding help
+
+Documentation can be found on [readthedocs](readthedocs.org/) under [micropython-ulab](https://micropython-ulab.readthedocs.io/en/latest), as well as at [circuitpython-ulab](
+https://circuitpython.readthedocs.io/en/latest/shared-bindings/ulab/__init__.html). A number of practical examples are listed in the excellent [circuitpython-ulab](https://learn.adafruit.com/ulab-crunch-numbers-fast-with-circuitpython/overview) overview.
 
 # Firmware
 
-Firmware for pyboard.v.1.1, and PYBD_SF6 is updated once in a while, and can be downloaded 
-from https://github.com/v923z/micropython-ulab/releases.
+Firmware for pyboard.v.1.1, and PYBD_SF6 is released once in a while, and can be downloaded
+from https://github.com/v923z/micropython-ulab/releases. Since a number of features can be
+set in the firmware (threading, support for SD card, LEDs, user switch etc.), and it is
+impossible to create something that suits everyone, these releases should only be used for
+quick testing of `ulab`. Otherwise, compilation from the source is required with
+the appropriate settings, which are usually defined in the `mpconfigboard.h` file of the port
+in question.
+
+`ulab` is also included in most development builds of `CircuitPython` for SAMD51 and nRF microcontrollers.
 
 ## Compiling
 
-If you want to try the latest version of `ulab`, or your hardware is 
+If you want to try the latest version of `ulab` on micropython, or your hardware is 
 different to pyboard.v.1.1, or PYBD_SF6, the firmware can be compiled 
 from the source by following these steps:
 
-First, you have to clone the micropython repository by running 
+First, you have to clone the `micropython` repository by running
 
 ```
 git clone https://github.com/micropython/micropython.git
@@ -30,18 +89,13 @@ on the command line. This will create a new repository with the name `micropytho
 git clone https://github.com/v923z/micropython-ulab.git ulab
 ```
 
-Then you have to include `ulab` in the compilation process by editing `mpconfigport.h` of the directory of the port for which you want to compile, so, still on the command line, navigate to `micropython/ports/unix`, or `micropython/ports/stm32`, or whichever port is your favourite, and edit the `mpconfigport.h` file there. All you have to do is add a single line at the end: 
-
-```
-#define MODULE_ULAB_ENABLED (1)
-```
-
-This line will inform the compiler that you want `ulab` in the resulting firmware. If you don't have the cross-compiler installed, your might want to do that now, for instance on Linux by executing 
+If you don't have the cross-compiler installed, your might want to do that now, for instance on Linux by executing
 
 ```
 sudo apt-get install gcc-arm-none-eabi
 ```
-If that was successful, you can try to run the make command in the port's directory as 
+
+If this step was successful, you can try to run the `make` command in the port's directory as
 ```
 make BOARD=PYBV11 USER_C_MODULES=../../../ulab all
 ```
@@ -49,8 +103,9 @@ which will prepare the firmware for pyboard.v.11. Similarly,
 ```
 make BOARD=PYBD_SF6 USER_C_MODULES=../../../ulab all
 ```
-will compile for the SF6 member of the PYBD series. Provided that you managed to compile the firmware, you would upload that by running
-either
+will compile for the SF6 member of the PYBD series. If your target is `unix`, you don't need to specify the `BOARD` parameter.
+
+Provided that you managed to compile the firmware, you would upload that by running either
 ```
 dfu-util --alt 0 -D firmware.dfu
 ```

build.sh

@@ -0,0 +1,17 @@
+#!/bin/sh
+set -e
+HERE="$(dirname -- "$(readlink -f -- "${0}")" )"
+[ -e micropython/py/py.mk ] || git clone https://github.com/micropython/micropython
+[ -e micropython/lib/libffi/autogen.sh ] || (cd micropython && git submodule update --init lib/libffi )
+#git clone https://github.com/micropython/micropython
+make -C micropython/mpy-cross -j$(nproc)
+make -C micropython/ports/unix -j$(nproc) deplibs
+make -C micropython/ports/unix -j$(nproc) USER_C_MODULES="${HERE}" DEBUG=1 STRIP=:
+
+if ! env MICROPY_MICROPYTHON=micropython/ports/unix/micropython micropython/tests/run-tests -d tests; then
+    for exp in *.exp; do
+        testbase=$(basename $exp .exp);
+        echo -e "\nFAILURE $testbase";
+        diff -u $testbase.exp $testbase.out;
+    done
+fi

code/approx.c

@@ -0,0 +1,396 @@
+/*
+ * This file is part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2020 Zoltán Vörös
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "py/obj.h"
+#include "py/runtime.h"
+#include "py/misc.h"
+#include "linalg.h"
+#include "approx.h"
+
+#if ULAB_APPROX_MODULE
+
+const mp_obj_float_t xtolerance = {{&mp_type_float}, 2.4e-7};
+const mp_obj_float_t rtolerance = {{&mp_type_float}, 0.0};
+
+STATIC mp_float_t approx_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(type->call(fun, nparams+1, 0, fargs));
+}
+
+STATIC mp_obj_t approx_bisect(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+	// Simple bisection routine
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} },
+        { MP_QSTR_xtol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} },
+        { MP_QSTR_maxiter, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 100} },
+    };
+
+    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(type->call == NULL) {
+		mp_raise_TypeError(translate("first argument must be a function"));
+	}
+	mp_float_t xtol = mp_obj_get_float(args[3].u_obj);
+	mp_obj_t *fargs = m_new(mp_obj_t, 1);
+	mp_float_t left, right;
+	mp_float_t x_mid;
+	mp_float_t a = mp_obj_get_float(args[1].u_obj);
+	mp_float_t b = mp_obj_get_float(args[2].u_obj);
+	left = approx_python_call(type, fun, a, fargs, 0);
+	right = approx_python_call(type, fun, b, fargs, 0);
+	if(left * right > 0) {
+		mp_raise_ValueError(translate("function has the same sign at the ends of interval"));
+	}
+	mp_float_t rtb = left < 0.0 ? a : b;
+	mp_float_t dx = left < 0.0 ? b - a : a - b;
+	
+	for(uint16_t i=0; i < args[4].u_int; i++) {
+		dx *= 0.5;
+		x_mid = rtb + dx;
+		if(approx_python_call(type, fun, x_mid, fargs, 0) < 0.0) {
+			rtb = x_mid;
+		}
+		if(MICROPY_FLOAT_C_FUN(fabs)(dx) < xtol) break;
+	}
+	return mp_obj_new_float(rtb);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(approx_bisect_obj, 3, approx_bisect);
+
+STATIC mp_obj_t approx_newton(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+	// this is actually the secant method, as the first derivative of the function 
+	// is not accepted as an argument. The function whose root we want to solve for 
+	// must depend on a single variable without parameters, i.e., f(x)
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_tol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} },
+        { MP_QSTR_rtol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&rtolerance)} },
+        { MP_QSTR_maxiter, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 50} },
+    };
+
+    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(type->call == NULL) {
+		mp_raise_TypeError(translate("first argument must be a function"));
+	}
+	mp_float_t x = mp_obj_get_float(args[1].u_obj);
+	mp_float_t tol = mp_obj_get_float(args[2].u_obj);
+	mp_float_t rtol = mp_obj_get_float(args[3].u_obj);
+	mp_float_t dx, df, fx;
+	dx = x > 0.0 ? APPROX_EPS * x : -APPROX_EPS * x;
+	mp_obj_t *fargs = m_new(mp_obj_t, 1);
+	for(uint16_t i=0; i < args[4].u_int; i++) {
+		fx = approx_python_call(type, fun, x, fargs, 0);
+		df = (approx_python_call(type, fun, x + dx, fargs, 0) - fx) / dx;
+		dx = fx / df;
+		x -= dx;
+		if(MICROPY_FLOAT_C_FUN(fabs)(dx) < (tol + rtol * MICROPY_FLOAT_C_FUN(fabs)(x))) break;
+	}
+	return mp_obj_new_float(x);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(approx_newton_obj, 2, approx_newton);
+
+STATIC mp_obj_t approx_fmin(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+	// downhill simplex method in 1D
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_xatol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} },
+        { MP_QSTR_fatol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} },
+        { MP_QSTR_maxiter, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} },
+    };
+
+    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(type->call == NULL) {
+		mp_raise_TypeError(translate("first argument must be a function"));
+	}
+	
+	// parameters controlling convergence conditions
+	mp_float_t xatol = mp_obj_get_float(args[2].u_obj);
+	mp_float_t fatol = mp_obj_get_float(args[3].u_obj);	
+	uint16_t maxiter = args[4].u_obj == mp_const_none ? 200 : mp_obj_get_int(args[4].u_obj);
+	
+	mp_float_t x0 = mp_obj_get_float(args[1].u_obj);
+	mp_float_t x1 = x0 != 0.0 ? (1.0 + APPROX_NONZDELTA) * x0 : APPROX_ZDELTA;
+	mp_obj_t *fargs = m_new(mp_obj_t, 1);
+	mp_float_t f0 = approx_python_call(type, fun, x0, fargs, 0);
+	mp_float_t f1 = approx_python_call(type, fun, x1, fargs, 0);
+	if(f1 < f0) {
+		SWAP(mp_float_t, x0, x1);
+		SWAP(mp_float_t, f0, f1);
+	}
+	for(uint16_t i=0; i < maxiter; i++) {
+		uint8_t shrink = 0;
+		f0 = approx_python_call(type, fun, x0, fargs, 0);
+		f1 = approx_python_call(type, fun, x1, fargs, 0);
+		
+		// reflection
+		mp_float_t xr = (1.0 + APPROX_ALPHA) * x0 - APPROX_ALPHA * x1;
+		mp_float_t fr = approx_python_call(type, fun, xr, fargs, 0);
+		if(fr < f0) { // expansion
+			mp_float_t xe = (1 + APPROX_ALPHA * APPROX_BETA) * x0 - APPROX_ALPHA * APPROX_BETA * x1;
+			mp_float_t fe = approx_python_call(type, fun, xe, fargs, 0);
+			if(fe < fr) {
+				x1 = xe;
+				f1 = fe;
+			} else {
+				x1 = xr;
+				f1 = fr;
+			}
+		} else {
+			if(fr < f1) { // contraction
+				mp_float_t xc = (1 + APPROX_GAMMA * APPROX_ALPHA) * x0 - APPROX_GAMMA * APPROX_ALPHA * x1;
+				mp_float_t fc = approx_python_call(type, fun, xc, fargs, 0);
+				if(fc < fr) {
+					x1 = xc;
+					f1 = fc;
+				} else {
+					shrink = 1;
+				}
+			} else { // inside contraction
+				mp_float_t xc = (1.0 - APPROX_GAMMA) * x0 + APPROX_GAMMA * x1;
+				mp_float_t fc = approx_python_call(type, fun, xc, fargs, 0);
+				if(fc < f1) {
+					x1 = xc;
+					f1 = fc;
+				} else {
+					shrink = 1;
+				}
+			}
+			if(shrink == 1) {
+				x1 = x0 + APPROX_DELTA * (x1 - x0);
+				f1 = approx_python_call(type, fun, x1, fargs, 0);
+			}
+			if((MICROPY_FLOAT_C_FUN(fabs)(f1 - f0) < fatol) || 
+				(MICROPY_FLOAT_C_FUN(fabs)(x1 - x0) < xatol)) {
+				break;
+			}
+			if(f1 < f0) {
+				SWAP(mp_float_t, x0, x1);
+				SWAP(mp_float_t, f0, f1);
+			}
+		}
+	}
+	return mp_obj_new_float(x0);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(approx_fmin_obj, 2, approx_fmin);
+
+#if 0
+static void approx_jacobi(const mp_obj_type_t *type, mp_obj_t fun, mp_float_t *x, mp_float_t *y, uint16_t len, mp_float_t *params, uint8_t nparams, mp_float_t *jacobi, mp_float_t *grad) {
+    /* Calculates the Jacobian and the gradient of the cost function
+     * 
+     * The entries in the Jacobian are 
+     * J(m, n) = de_m/da_n, 
+     * 
+     * where 
+     * 
+     * e_m = (f(x_m, a1, a2, ...) - y_m)/sigma_m is the error at x_m, 
+     * 
+     * and
+     * 
+     * a1, a2, ..., a_n are the free parameters
+     */
+    mp_obj_t *fargs0 = m_new(mp_obj_t, lenp+1); 
+    mp_obj_t *fargs1 = m_new(mp_obj_t, lenp+1);
+    for(uint8_t p=0; p < nparams; p++) {
+        fargs0[p+1] = mp_obj_new_float(params[p]);
+        fargs1[p+1] = mp_obj_new_float(params[p]);
+    }
+    for(uint8_t p=0; p < nparams; p++) {
+        mp_float_t da = params[p] != 0.0 ? (1.0 + APPROX_NONZDELTA) * params[p] : APPROX_ZDELTA;
+        fargs1[p+1] = mp_obj_new_float(params[p] + da);
+        grad[p] = 0.0;
+        for(uint16_t i=0; i < len; i++) {
+            mp_float_t f0 = approx_python_call(type, fun, x[i], fargs0, nparams);
+            mp_float_t f1 = approx_python_call(type, fun, x[i], fargs1, nparams);
+            jacobi[i*nparamp+p] = (f1 - f0) / da;
+            grad[p] += (f0 - y[i]) * jacobi[i*nparamp+p];
+        }
+        fargs1[p+1] = fargs0[p+1]; // set back to the original value
+    }
+}
+
+static void approx_delta(mp_float_t *jacobi, mp_float_t *grad, uint16_t len, uint8_t nparams, mp_float_t lambda) {
+    // 
+}
+
+mp_obj_t approx_curve_fit(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+	// Levenberg-Marquardt non-linear fit
+    // The implementation follows the introductory discussion in Mark Tanstrum's paper, https://arxiv.org/abs/1201.5885
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_p0, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_xatol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} },
+        { MP_QSTR_fatol, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_PTR(&xtolerance)} },
+        { MP_QSTR_maxiter, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} },
+    };
+
+    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(type->call == NULL) {
+		mp_raise_TypeError(translate("first argument must be a function"));
+	}
+    
+	mp_obj_t x_obj = args[1].u_obj;
+	mp_obj_t y_obj = args[2].u_obj;
+    mp_obj_t p0_obj = args[3].u_obj;
+    if(!ndarray_object_is_nditerable(x_obj) || !ndarray_object_is_nditerable(y_obj)) {
+        mp_raise_TypeError(translate("data must be iterable"));
+    }
+    if(!ndarray_object_is_nditerable(p0_obj)) {
+        mp_raise_TypeError(translate("initial values must be iterable"));
+    }
+    uint16_t len = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(x_obj));
+    uint8_t lenp = (uint8_t)mp_obj_get_int(mp_obj_len_maybe(p0_obj));
+    if(len != (uint16_t)mp_obj_get_int(mp_obj_len_maybe(y_obj))) {
+        mp_raise_ValueError(translate("data must be of equal length"));
+    }
+
+    mp_float_t *x = m_new(mp_float_t, len);
+    fill_array_iterable(x, x_obj);
+    mp_float_t *y = m_new(mp_float_t, len);
+    fill_array_iterable(y, y_obj);
+    mp_float_t *p0 = m_new(mp_float_t, lenp);
+    fill_array_iterable(p0, p0_obj);
+    mp_float_t *grad = m_new(mp_float_t, len);
+    mp_float_t *jacobi = m_new(mp_float_t, len*len);
+    mp_obj_t *fargs = m_new(mp_obj_t, lenp+1);
+
+    m_del(mp_float_t, p0, lenp);
+	// parameters controlling convergence conditions
+	//mp_float_t xatol = mp_obj_get_float(args[2].u_obj);
+	//mp_float_t fatol = mp_obj_get_float(args[3].u_obj);	
+    
+    // this has finite binary representation; we will multiply/divide by 4
+    //mp_float_t lambda = 0.0078125; 
+    
+    //linalg_invert_matrix(mp_float_t *data, size_t N)
+    
+    m_del(mp_float_t, x, len);
+    m_del(mp_float_t, y, len);
+    m_del(mp_float_t, grad, len);
+    m_del(mp_float_t, jacobi, len*len);
+    m_del(mp_obj_t, fargs, lenp+1);
+    return mp_const_none;
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(approx_curve_fit_obj, 2, approx_curve_fit);
+
+#endif
+
+STATIC mp_obj_t approx_interp(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_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_left, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} },
+        { MP_QSTR_right, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} },
+    };
+	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);
+	
+	ndarray_obj_t *x = ndarray_from_mp_obj(args[0].u_obj);
+	ndarray_obj_t *xp = ndarray_from_mp_obj(args[1].u_obj); // xp must hold an increasing sequence of independent values
+	ndarray_obj_t *fp = ndarray_from_mp_obj(args[2].u_obj);
+	// TODO: check if the shape is (1, n), or (m, 1)
+	if(((xp->m != 1) && (xp->n != 1)) || ((fp->m != 1) && (fp->n != 1)) || 
+		(xp->array->len < 2) || (fp->array->len < 2) || (xp->array->len != fp->array->len)) {
+		mp_raise_ValueError(translate("interp is defined for 1D arrays of equal length"));
+	}
+	ndarray_obj_t *y = create_new_ndarray(x->m, x->n, NDARRAY_FLOAT);
+	mp_float_t left_value, right_value;
+	mp_float_t xp_left = ndarray_get_float_value(xp->array->items, xp->array->typecode, 0);
+	mp_float_t xp_right = ndarray_get_float_value(xp->array->items, xp->array->typecode, xp->array->len-1);
+	if(args[3].u_obj == mp_const_none) {
+		left_value = ndarray_get_float_value(fp->array->items, fp->array->typecode, 0);
+	} else {
+		left_value = mp_obj_get_float(args[3].u_obj);
+	}
+	if(args[4].u_obj == mp_const_none) {
+		right_value = ndarray_get_float_value(fp->array->items, fp->array->typecode, fp->array->len-1);
+	} else {
+		right_value = mp_obj_get_float(args[4].u_obj);
+	}
+	mp_float_t *yarray = (mp_float_t *)y->array->items;
+	for(size_t i=0; i < x->array->len; i++, yarray++) {
+		mp_float_t x_value = ndarray_get_float_value(x->array->items, x->array->typecode, i);
+		if(x_value <= xp_left) {
+			*yarray = left_value;
+		} else if(x_value >= xp_right) {
+			*yarray = right_value;
+		} else { // do the binary search here
+			mp_float_t xp_left_, xp_right_;
+			mp_float_t fp_left, fp_right;
+			size_t left_index = 0, right_index = xp->array->len - 1, middle_index;
+			while(right_index - left_index > 1) {
+				middle_index = left_index + (right_index - left_index) / 2;
+				mp_float_t xp_middle = ndarray_get_float_value(xp->array->items, xp->array->typecode, middle_index);
+				if(x_value <= xp_middle) {
+					right_index = middle_index;
+				} else {
+					left_index = middle_index;
+				}
+			}
+			xp_left_ = ndarray_get_float_value(xp->array->items, xp->array->typecode, left_index);
+			xp_right_ = ndarray_get_float_value(xp->array->items, xp->array->typecode, right_index);
+			fp_left = ndarray_get_float_value(fp->array->items, fp->array->typecode, left_index);
+			fp_right = ndarray_get_float_value(fp->array->items, fp->array->typecode, right_index);
+			*yarray = fp_left + (x_value - xp_left_) * (fp_right - fp_left) / (xp_right_ - xp_left_);
+		}
+	}
+	return MP_OBJ_FROM_PTR(y);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(approx_interp_obj, 2, approx_interp);
+
+STATIC const mp_rom_map_elem_t ulab_approx_globals_table[] = {
+    { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_approx) },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_bisect), (mp_obj_t)&approx_bisect_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_newton), (mp_obj_t)&approx_newton_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_fmin), (mp_obj_t)&approx_fmin_obj },
+//    { MP_OBJ_NEW_QSTR(MP_QSTR_curve_fit), (mp_obj_t)&approx_curve_fit_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_interp), (mp_obj_t)&approx_interp_obj },
+};
+
+STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_approx_globals, ulab_approx_globals_table);
+
+mp_obj_module_t ulab_approx_module = {
+    .base = { &mp_type_module },
+    .globals = (mp_obj_dict_t*)&mp_module_ulab_approx_globals,
+};
+
+#endif
+

code/approx.h

@@ -0,0 +1,42 @@
+
+/*
+ * This file is part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2020 Zoltán Vörös
+*/
+
+#ifndef _APPROX_
+#define _APPROX_
+
+#include "ulab.h"
+#include "ndarray.h"
+
+#if ULAB_APPROX_MODULE
+
+// this typedef is lifted from objfloat.c, because mp_obj_float_t is not exposed
+typedef struct _mp_obj_float_t {
+    mp_obj_base_t base;
+    mp_float_t value;
+} mp_obj_float_t;
+
+#define 	APPROX_EPS			1.0e-4
+#define 	APPROX_NONZDELTA	0.05
+#define 	APPROX_ZDELTA		0.00025
+#define		APPROX_ALPHA		1.0
+#define		APPROX_BETA			2.0
+#define		APPROX_GAMMA		0.5
+#define		APPROX_DELTA		0.5
+
+extern mp_obj_module_t ulab_approx_module;
+
+MP_DECLARE_CONST_FUN_OBJ_KW(approx_bisect_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(approx_newton_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(approx_fmin_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(approx_interp_obj);
+
+#endif
+#endif

code/compare.c

@@ -0,0 +1,194 @@
+
+/*
+ * This file is part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2020 Zoltán Vörös
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "py/obj.h"
+#include "py/runtime.h"
+#include "py/misc.h"
+#include "compare.h"
+
+#if ULAB_COMPARE_MODULE
+
+static mp_obj_t compare_function(mp_obj_t x1, mp_obj_t x2, uint8_t comptype) {
+	// the function is implemented for scalars and ndarrays only, with partial 
+	// broadcasting: arguments must be either scalars, or ndarrays of equal size/shape
+	if(!(MP_OBJ_IS_INT(x1) || mp_obj_is_float(x1) || MP_OBJ_IS_TYPE(x1, &ulab_ndarray_type)) &&
+		!(MP_OBJ_IS_INT(x2) || mp_obj_is_float(x2) || MP_OBJ_IS_TYPE(x2, &ulab_ndarray_type))) {
+		mp_raise_TypeError(translate("function is implemented for scalars and ndarrays only"));
+	}
+	ndarray_obj_t *ndarray1 = ndarray_from_mp_obj(x1);
+	ndarray_obj_t *ndarray2 = ndarray_from_mp_obj(x2); 
+	// check, whether partial broadcasting is possible here
+	if((ndarray1->m != ndarray2->m) || (ndarray1->n != ndarray2->n)) {
+		if((ndarray1->array->len != 1) && (ndarray2->array->len != 1)) {
+            mp_raise_ValueError(translate("operands could not be broadcast together"));
+		}
+	}
+	if((comptype == MP_BINARY_OP_EQUAL) || (comptype == MP_BINARY_OP_NOT_EQUAL)) {
+		return ndarray_binary_op(comptype, x1, x2);
+	}
+	size_t m = MAX(ndarray1->m, ndarray2->m);
+	size_t n = MAX(ndarray1->n, ndarray2->n);
+	size_t len = MAX(ndarray1->array->len, ndarray2->array->len);
+	size_t inc1 = ndarray1->array->len == 1 ? 0 : 1;
+	size_t inc2 = ndarray2->array->len == 1 ? 0 : 1;
+	// These are the upcasting rules
+	// float always becomes float
+	// operation on identical types preserves type
+	// uint8 + int8 => int16
+	// uint8 + int16 => int16
+	// uint8 + uint16 => uint16
+	// int8 + int16 => int16
+	// int8 + uint16 => uint16
+	// uint16 + int16 => float
+	// The parameters of RUN_BINARY_LOOP are 
+	// typecode of result, type_out, type_left, type_right, lhs operand, rhs operand, operator
+	if(ndarray1->array->typecode == NDARRAY_UINT8) {
+		if(ndarray2->array->typecode == NDARRAY_UINT8) {
+			RUN_COMPARE_LOOP(NDARRAY_UINT8, uint8_t, uint8_t, uint8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT8) {
+			RUN_COMPARE_LOOP(NDARRAY_INT16, int16_t, uint8_t, int8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_UINT16) {
+			RUN_COMPARE_LOOP(NDARRAY_UINT16, uint16_t, uint8_t, uint16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT16) {
+			RUN_COMPARE_LOOP(NDARRAY_INT16, int16_t, uint8_t, int16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_FLOAT) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		}
+	} else if(ndarray1->array->typecode == NDARRAY_INT8) {
+		if(ndarray2->array->typecode == NDARRAY_UINT8) {
+			RUN_COMPARE_LOOP(NDARRAY_INT16, int16_t, int8_t, uint8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT8) {
+			RUN_COMPARE_LOOP(NDARRAY_INT8, int8_t, int8_t, int8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_UINT16) {
+			RUN_COMPARE_LOOP(NDARRAY_INT16, int16_t, int8_t, uint16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT16) {
+			RUN_COMPARE_LOOP(NDARRAY_INT16, int16_t, int8_t, int16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_FLOAT) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, int8_t, mp_float_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		}                
+	} else if(ndarray1->array->typecode == NDARRAY_UINT16) {
+		if(ndarray2->array->typecode == NDARRAY_UINT8) {
+			RUN_COMPARE_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT8) {
+			RUN_COMPARE_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, int8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_UINT16) {
+			RUN_COMPARE_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT16) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, int16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_FLOAT) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		}
+	} else if(ndarray1->array->typecode == NDARRAY_INT16) {
+		if(ndarray2->array->typecode == NDARRAY_UINT8) {
+			RUN_COMPARE_LOOP(NDARRAY_INT16, int16_t, int16_t, uint8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT8) {
+			RUN_COMPARE_LOOP(NDARRAY_INT16, int16_t, int16_t, int8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_UINT16) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, int16_t, uint16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT16) {
+			RUN_COMPARE_LOOP(NDARRAY_INT16, int16_t, int16_t, int16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_FLOAT) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, mp_float_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		}
+	} else if(ndarray1->array->typecode == NDARRAY_FLOAT) {
+		if(ndarray2->array->typecode == NDARRAY_UINT8) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT8) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int8_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_UINT16) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_INT16) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int16_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		} else if(ndarray2->array->typecode == NDARRAY_FLOAT) {
+			RUN_COMPARE_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, mp_float_t, ndarray1, ndarray2, comptype, m, n, len, inc1, inc2);
+		}
+	}
+    return mp_const_none; // we should never reach this point
+}
+
+static mp_obj_t compare_equal_helper(mp_obj_t x1, mp_obj_t x2, uint8_t comptype) {
+	// scalar comparisons should return a single object of mp_obj_t type
+	mp_obj_t result = compare_function(x1, x2, comptype);
+	if((MP_OBJ_IS_INT(x1) || mp_obj_is_float(x1)) && (MP_OBJ_IS_INT(x2) || mp_obj_is_float(x2))) {
+		mp_obj_iter_buf_t iter_buf;
+		mp_obj_t iterable = mp_getiter(result, &iter_buf);
+		mp_obj_t item = mp_iternext(iterable);
+		return item;
+	}
+	return result;	
+
+}
+static mp_obj_t compare_equal(mp_obj_t x1, mp_obj_t x2) {
+	return compare_equal_helper(x1, x2, MP_BINARY_OP_EQUAL);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_2(compare_equal_obj, compare_equal);
+
+static mp_obj_t compare_not_equal(mp_obj_t x1, mp_obj_t x2) {
+	return compare_equal_helper(x1, x2, MP_BINARY_OP_NOT_EQUAL);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_2(compare_not_equal_obj, compare_not_equal);
+
+static mp_obj_t compare_minimum(mp_obj_t x1, mp_obj_t x2) {
+	// extra round, so that we can return minimum(3, 4) properly
+	mp_obj_t result = compare_function(x1, x2, COMPARE_MINIMUM);
+	if((MP_OBJ_IS_INT(x1) || mp_obj_is_float(x1)) && (MP_OBJ_IS_INT(x2) || mp_obj_is_float(x2))) {
+		ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(result);
+		return mp_binary_get_val_array(ndarray->array->typecode, ndarray->array->items, 0);
+	}
+	return result;
+}
+
+MP_DEFINE_CONST_FUN_OBJ_2(compare_minimum_obj, compare_minimum);
+
+static mp_obj_t compare_maximum(mp_obj_t x1, mp_obj_t x2) {
+	// extra round, so that we can return maximum(3, 4) properly
+	mp_obj_t result = compare_function(x1, x2, COMPARE_MAXIMUM);
+	if((MP_OBJ_IS_INT(x1) || mp_obj_is_float(x1)) && (MP_OBJ_IS_INT(x2) || mp_obj_is_float(x2))) {
+		ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(result);
+		return mp_binary_get_val_array(ndarray->array->typecode, ndarray->array->items, 0);
+	}
+	return result;
+}
+
+MP_DEFINE_CONST_FUN_OBJ_2(compare_maximum_obj, compare_maximum);
+
+static mp_obj_t compare_clip(mp_obj_t x1, mp_obj_t x2, mp_obj_t x3) {
+	// Note: this function could be made faster by implementing a single-loop comparison in 
+	// RUN_COMPARE_LOOP. However, that would add around 2 kB of compile size, while we 
+	// would not gain a factor of two in speed, since the two comparisons should still be 
+	// evaluated. In contrast, calling the function twice adds only 140 bytes to the firmware
+	return compare_function(x2, compare_function(x1, x3, COMPARE_MINIMUM), COMPARE_MAXIMUM);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_3(compare_clip_obj, compare_clip);
+
+STATIC const mp_rom_map_elem_t ulab_compare_globals_table[] = {
+    { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_compare) },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_equal), (mp_obj_t)&compare_equal_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_not_equal), (mp_obj_t)&compare_not_equal_obj },
+	{ MP_OBJ_NEW_QSTR(MP_QSTR_maximum), (mp_obj_t)&compare_maximum_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_minimum), (mp_obj_t)&compare_minimum_obj },
+	{ MP_OBJ_NEW_QSTR(MP_QSTR_clip), (mp_obj_t)&compare_clip_obj },
+};
+
+STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_compare_globals, ulab_compare_globals_table);
+
+mp_obj_module_t ulab_compare_module = {
+    .base = { &mp_type_module },
+    .globals = (mp_obj_dict_t*)&mp_module_ulab_compare_globals,
+};
+
+#endif

code/compare.h

@@ -0,0 +1,45 @@
+
+/*
+ * This file is part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2020 Zoltán Vörös
+*/
+
+#ifndef _COMPARE_
+#define _COMPARE_
+
+#include "ulab.h"
+#include "ndarray.h"
+
+#if ULAB_COMPARE_MODULE
+
+enum COMPARE_FUNCTION_TYPE {
+    COMPARE_MINIMUM,
+    COMPARE_MAXIMUM,
+    COMPARE_CLIP,
+};
+
+extern mp_obj_module_t ulab_compare_module;
+
+MP_DECLARE_CONST_FUN_OBJ_2(compare_equal_obj);
+MP_DECLARE_CONST_FUN_OBJ_2(compare_not_equal_obj);
+MP_DECLARE_CONST_FUN_OBJ_2(compare_minimum_obj);
+MP_DECLARE_CONST_FUN_OBJ_2(compare_maximum_obj);
+MP_DECLARE_CONST_FUN_OBJ_3(compare_clip_obj);
+
+#define RUN_COMPARE_LOOP(typecode, type_out, type1, type2, nd1, nd2, op, m, n, len, inc1, inc2) do {\
+	type1 *array1 = (type1 *)(nd1)->array->items;\
+    type2 *array2 = (type2 *)(nd2)->array->items;\
+	ndarray_obj_t *out = create_new_ndarray((m), (n), (typecode));\
+	type_out *(odata) = (type_out *)out->array->items;\
+	if((op) == COMPARE_MINIMUM) { for(size_t i=0; i < (len); i++, array1+=(inc1), array2+=(inc2)) *odata++ = *array1 < *array2 ? *array1 : *array2; }\
+	if((op) == COMPARE_MAXIMUM) { for(size_t i=0; i < (len); i++, array1+=(inc1), array2+=(inc2)) *odata++ = *array1 > *array2 ? *array1 : *array2; }\
+	return MP_OBJ_FROM_PTR(out);\
+} while(0)
+
+#endif
+#endif

code/create.c

@@ -0,0 +1,155 @@
+/*
+ * This file is part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2020 Jeff Epler for Adafruit Industries
+ * 				 2019-2020 Zoltán Vörös
+*/
+
+
+#include "py/obj.h"
+#include "py/runtime.h"
+#include "create.h"
+
+static mp_obj_t create_zeros_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t kind) {
+    static const mp_arg_t allowed_args[] = {
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_obj = MP_OBJ_NULL} } ,
+        { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },
+    };
+
+    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);
+    
+    uint8_t dtype = args[1].u_int;
+    if(!MP_OBJ_IS_INT(args[0].u_obj) && !MP_OBJ_IS_TYPE(args[0].u_obj, &mp_type_tuple)) {
+        mp_raise_TypeError(translate("input argument must be an integer or a 2-tuple"));
+    }
+    ndarray_obj_t *ndarray = NULL;
+    if(MP_OBJ_IS_INT(args[0].u_obj)) {
+        size_t n = mp_obj_get_int(args[0].u_obj);
+        ndarray = create_new_ndarray(1, n, dtype);
+    } else if(MP_OBJ_IS_TYPE(args[0].u_obj, &mp_type_tuple)) {
+        mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(args[0].u_obj);
+        if(tuple->len != 2) {
+            mp_raise_TypeError(translate("input argument must be an integer or a 2-tuple"));
+        }
+        ndarray = create_new_ndarray(mp_obj_get_int(tuple->items[0]), 
+                                                  mp_obj_get_int(tuple->items[1]), dtype);
+    }
+    if(kind == 1) {
+        mp_obj_t one = mp_obj_new_int(1);
+        for(size_t i=0; i < ndarray->array->len; i++) {
+            mp_binary_set_val_array(dtype, ndarray->array->items, i, one);
+        }
+    }
+    return MP_OBJ_FROM_PTR(ndarray);
+}
+
+mp_obj_t create_zeros(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+    return create_zeros_ones(n_args, pos_args, kw_args, 0);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(create_zeros_obj, 0, create_zeros);
+
+mp_obj_t create_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+    return create_zeros_ones(n_args, pos_args, kw_args, 1);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(create_ones_obj, 0, create_ones);
+
+mp_obj_t create_eye(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_INT, {.u_int = 0} },
+        { MP_QSTR_M, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_k, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} },        
+        { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },
+    };
+
+    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);
+
+    size_t n = args[0].u_int, m;
+    size_t k = args[2].u_int;
+    uint8_t dtype = args[3].u_int;
+    if(args[1].u_rom_obj == mp_const_none) {
+        m = n;
+    } else {
+        m = mp_obj_get_int(args[1].u_rom_obj);
+    }
+    
+    ndarray_obj_t *ndarray = create_new_ndarray(m, n, dtype);
+    mp_obj_t one = mp_obj_new_int(1);
+    size_t i = 0;
+    if((k >= 0) && (k < n)) {
+        while(k < n) {
+            mp_binary_set_val_array(dtype, ndarray->array->items, i*n+k, one);
+            k++;
+            i++;
+        }
+    } else if((k < 0) && (-k < m)) {
+        k = -k;
+        i = 0;
+        while(k < m) {
+            mp_binary_set_val_array(dtype, ndarray->array->items, k*n+i, one);
+            k++;
+            i++;
+        }
+    }
+    return MP_OBJ_FROM_PTR(ndarray);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(create_eye_obj, 0, create_eye);
+
+mp_obj_t create_linspace(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_const_none } },
+        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
+        { MP_QSTR_num, MP_ARG_INT, {.u_int = 50} },
+        { MP_QSTR_endpoint, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_true} },
+        { MP_QSTR_retstep, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_false} },
+        { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },
+    };
+
+    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);
+
+    uint16_t len = args[2].u_int;
+    if(len < 2) {
+        mp_raise_ValueError(translate("number of points must be at least 2"));
+    }
+    mp_float_t value, step;
+    value = mp_obj_get_float(args[0].u_obj);
+    uint8_t typecode = args[5].u_int;
+    if(args[3].u_obj == mp_const_true) step = (mp_obj_get_float(args[1].u_obj)-value)/(len-1);
+    else step = (mp_obj_get_float(args[1].u_obj)-value)/len;
+    ndarray_obj_t *ndarray = create_new_ndarray(1, len, typecode);
+    if(typecode == NDARRAY_UINT8) {
+        uint8_t *array = (uint8_t *)ndarray->array->items;
+        for(size_t i=0; i < len; i++, value += step) array[i] = (uint8_t)value;
+    } else if(typecode == NDARRAY_INT8) {
+        int8_t *array = (int8_t *)ndarray->array->items;
+        for(size_t i=0; i < len; i++, value += step) array[i] = (int8_t)value;
+    } else if(typecode == NDARRAY_UINT16) {
+        uint16_t *array = (uint16_t *)ndarray->array->items;
+        for(size_t i=0; i < len; i++, value += step) array[i] = (uint16_t)value;
+    } else if(typecode == NDARRAY_INT16) {
+        int16_t *array = (int16_t *)ndarray->array->items;
+        for(size_t i=0; i < len; i++, value += step) array[i] = (int16_t)value;
+    } else {
+        mp_float_t *array = (mp_float_t *)ndarray->array->items;
+        for(size_t i=0; i < len; i++, value += step) array[i] = value;
+    }
+    if(args[4].u_obj == mp_const_false) {
+        return MP_OBJ_FROM_PTR(ndarray);
+    } else {
+        mp_obj_t tuple[2];
+        tuple[0] = ndarray;
+        tuple[1] = mp_obj_new_float(step);
+        return mp_obj_new_tuple(2, tuple);
+    }
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(create_linspace_obj, 2, create_linspace);

code/create.h

@@ -0,0 +1,28 @@
+/*
+ * This file is part of the micropython-ulab project,
+ *
+ * https://github.com/v923z/micropython-ulab
+ *
+ * The MIT License (MIT)
+ *
+ * Copyright (c) 2020 Jeff Epler for Adafruit Industries
+ * 				 2019-2020 Zoltán Vörös
+*/
+
+#ifndef _CREATE_
+#define _CREATE_
+
+#include "ulab.h"
+#include "ndarray.h"
+
+mp_obj_t create_zeros(size_t , const mp_obj_t *, mp_map_t *);
+mp_obj_t create_ones(size_t , const mp_obj_t *, mp_map_t *);
+mp_obj_t create_eye(size_t , const mp_obj_t *, mp_map_t *);
+mp_obj_t create_linspace(size_t , const mp_obj_t *, mp_map_t *);
+
+MP_DECLARE_CONST_FUN_OBJ_KW(create_ones_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(create_zeros_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(create_eye_obj);
+MP_DECLARE_CONST_FUN_OBJ_KW(create_linspace_obj);
+
+#endif

code/extras.c

@@ -19,13 +19,24 @@
 
 #if ULAB_EXTRAS_MODULE
 
-STATIC const mp_rom_map_elem_t ulab_filter_globals_table[] = {
+static mp_obj_t extras_spectrogram(size_t n_args, const mp_obj_t *args) {
+    if(n_args == 2) {
+        return fft_fft_ifft_spectrum(n_args, args[0], args[1], FFT_SPECTRUM);
+    } else {
+        return fft_fft_ifft_spectrum(n_args, args[0], mp_const_none, FFT_SPECTRUM);
+    }
+}
+
+MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(extras_spectrogram_obj, 1, 2, extras_spectrogram);
+
+STATIC const mp_rom_map_elem_t ulab_extras_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_extras) },
+	{ MP_OBJ_NEW_QSTR(MP_QSTR_spectrogram), (mp_obj_t)&extras_spectrogram_obj },
 };
 
 STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_extras_globals, ulab_extras_globals_table);
 
-mp_obj_module_t ulab_filter_module = {
+mp_obj_module_t ulab_extras_module = {
     .base = { &mp_type_module },
     .globals = (mp_obj_dict_t*)&mp_module_ulab_extras_globals,
 };

code/extras.h

@@ -14,6 +14,7 @@
 
 #include "ulab.h"
 #include "ndarray.h"
+#include "fft.h"
 
 #if ULAB_EXTRAS_MODULE
 

code/fft.c

@@ -23,13 +23,7 @@
 
 #if ULAB_FFT_MODULE
 
-enum FFT_TYPE {
-    FFT_FFT,
-    FFT_IFFT,
-    FFT_SPECTRUM,
-};
-
-void fft_kernel(mp_float_t *real, mp_float_t *imag, int n, int isign) {
+static void fft_kernel(mp_float_t *real, mp_float_t *imag, int n, int isign) {
     // This is basically a modification of four1 from Numerical Recipes
     // The main difference is that this function takes two arrays, one 
     // for the real, and one for the imaginary parts. 
@@ -152,7 +146,7 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
     }
 }
 
-mp_obj_t fft_fft(size_t n_args, const mp_obj_t *args) {
+static mp_obj_t fft_fft(size_t n_args, const mp_obj_t *args) {
     if(n_args == 2) {
         return fft_fft_ifft_spectrum(n_args, args[0], args[1], FFT_FFT);
     } else {
@@ -162,7 +156,7 @@ mp_obj_t fft_fft(size_t n_args, const mp_obj_t *args) {
 
 MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj, 1, 2, fft_fft);
 
-mp_obj_t fft_ifft(size_t n_args, const mp_obj_t *args) {
+static mp_obj_t fft_ifft(size_t n_args, const mp_obj_t *args) {
     if(n_args == 2) {
         return fft_fft_ifft_spectrum(n_args, args[0], args[1], FFT_IFFT);
     } else {
@@ -172,22 +166,10 @@ mp_obj_t fft_ifft(size_t n_args, const mp_obj_t *args) {
 
 MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_ifft_obj, 1, 2, fft_ifft);
 
-mp_obj_t fft_spectrum(size_t n_args, const mp_obj_t *args) {
-    if(n_args == 2) {
-        return fft_fft_ifft_spectrum(n_args, args[0], args[1], FFT_SPECTRUM);
-    } else {
-        return fft_fft_ifft_spectrum(n_args, args[0], mp_const_none, FFT_SPECTRUM);
-    }
-}
-
-MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_spectrum_obj, 1, 2, fft_spectrum);
-
-#if !CIRCUITPY
 STATIC const mp_rom_map_elem_t ulab_fft_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_fft) },
     { MP_OBJ_NEW_QSTR(MP_QSTR_fft), (mp_obj_t)&fft_fft_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_ifft), (mp_obj_t)&fft_ifft_obj },
-    { MP_OBJ_NEW_QSTR(MP_QSTR_spectrum), (mp_obj_t)&fft_spectrum_obj },
 };
 
 STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_fft_globals, ulab_fft_globals_table);
@@ -196,6 +178,5 @@ mp_obj_module_t ulab_fft_module = {
     .base = { &mp_type_module },
     .globals = (mp_obj_dict_t*)&mp_module_ulab_fft_globals,
 };
-#endif
 
 #endif

code/fft.h

@@ -19,6 +19,12 @@
 
 #define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }
 
+enum FFT_TYPE {
+    FFT_FFT,
+    FFT_IFFT,
+    FFT_SPECTRUM,
+};
+
 #if ULAB_FFT_MODULE
 
 extern mp_obj_module_t ulab_fft_module;
@@ -27,5 +33,7 @@ MP_DECLARE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj);
 MP_DECLARE_CONST_FUN_OBJ_VAR_BETWEEN(fft_ifft_obj);
 MP_DECLARE_CONST_FUN_OBJ_VAR_BETWEEN(fft_spectrum_obj);
 
+mp_obj_t fft_fft_ifft_spectrum(size_t , mp_obj_t , mp_obj_t , uint8_t );
+
 #endif
 #endif

code/filter.c

@@ -18,14 +18,14 @@
 #include "filter.h"
 
 #if ULAB_FILTER_MODULE
-mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     static const mp_arg_t allowed_args[] = {
         { MP_QSTR_a, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
         { MP_QSTR_v, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
 
     if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type) || !MP_OBJ_IS_TYPE(args[1].u_obj, &ulab_ndarray_type)) {
         mp_raise_TypeError(translate("convolve arguments must be ndarrays"));
@@ -33,8 +33,8 @@ mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a
 
     ndarray_obj_t *a = MP_OBJ_TO_PTR(args[0].u_obj);
     ndarray_obj_t *c = MP_OBJ_TO_PTR(args[1].u_obj);
-    int len_a = a->array->len;
-    int len_c = c->array->len;
+    size_t len_a = a->array->len;
+    size_t len_c = c->array->len;
     // deal with linear arrays only
     if(a->m*a->n != len_a || c->m*c->n != len_c) {
         mp_raise_TypeError(translate("convolve arguments must be linear arrays"));
@@ -84,7 +84,6 @@ mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a
 
 MP_DEFINE_CONST_FUN_OBJ_KW(filter_convolve_obj, 2, filter_convolve);
 
-#if !CIRCUITPY
 STATIC const mp_rom_map_elem_t ulab_filter_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_filter) },
     { MP_OBJ_NEW_QSTR(MP_QSTR_convolve), (mp_obj_t)&filter_convolve_obj },
@@ -96,6 +95,5 @@ mp_obj_module_t ulab_filter_module = {
     .base = { &mp_type_module },
     .globals = (mp_obj_dict_t*)&mp_module_ulab_filter_globals,
 };
-#endif
 
 #endif

code/linalg.c

@@ -19,14 +19,27 @@
 
 #if ULAB_LINALG_MODULE
 
-mp_obj_t linalg_size(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static ndarray_obj_t *linalg_object_is_square(mp_obj_t obj) {
+	// Returns an ndarray, if the object is a square ndarray, 
+	// raises the appropriate exception otherwise
+	if(!MP_OBJ_IS_TYPE(obj, &ulab_ndarray_type)) {
+        mp_raise_TypeError(translate("size is defined for ndarrays only"));
+    }
+	ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(obj);
+	if(ndarray->m != ndarray->n) {
+		mp_raise_ValueError(translate("input must be square matrix"));
+	}
+	return ndarray;
+}
+
+static mp_obj_t linalg_size(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_const_none } },
         { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
 
     if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {
         mp_raise_TypeError(translate("size is defined for ndarrays only"));
@@ -101,18 +114,8 @@ bool linalg_invert_matrix(mp_float_t *data, size_t N) {
     return true;
 }
 
-mp_obj_t linalg_inv(mp_obj_t o_in) {
-    // since inv is not a class method, we have to inspect the input argument first
-    if(!MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {
-        mp_raise_TypeError(translate("only ndarrays can be inverted"));
-    }
-    ndarray_obj_t *o = MP_OBJ_TO_PTR(o_in);
-    if(!MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {
-        mp_raise_TypeError(translate("only ndarray objects can be inverted"));
-    }
-    if(o->m != o->n) {
-        mp_raise_ValueError(translate("only square matrices can be inverted"));
-    }
+static mp_obj_t linalg_inv(mp_obj_t o_in) {
+	ndarray_obj_t *o = linalg_object_is_square(o_in);
     ndarray_obj_t *inverted = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT);
     mp_float_t *data = (mp_float_t *)inverted->array->items;
     mp_obj_t elem;
@@ -136,7 +139,7 @@ mp_obj_t linalg_inv(mp_obj_t o_in) {
 
 MP_DEFINE_CONST_FUN_OBJ_1(linalg_inv_obj, linalg_inv);
 
-mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {
+static mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {
     // TODO: should the results be upcast?
     if(!MP_OBJ_IS_TYPE(_m1, &ulab_ndarray_type) || !MP_OBJ_IS_TYPE(_m2, &ulab_ndarray_type)) {
         mp_raise_TypeError(translate("arguments must be ndarrays"));
@@ -149,17 +152,16 @@ mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {
     // TODO: numpy uses upcasting here
     ndarray_obj_t *out = create_new_ndarray(m1->m, m2->n, NDARRAY_FLOAT);
     mp_float_t *outdata = (mp_float_t *)out->array->items;
-    mp_float_t sum, v1, v2;
     for(size_t i=0; i < m1->m; i++) { // rows of m1
         for(size_t j=0; j < m2->n; j++) { // columns of m2
-            sum = 0.0;
+            mp_float_t sum = 0.0, v1, v2;
             for(size_t k=0; k < m2->m; k++) {
                 // (i, k) * (k, j)
                 v1 = ndarray_get_float_value(m1->array->items, m1->array->typecode, i*m1->n+k);
                 v2 = ndarray_get_float_value(m2->array->items, m2->array->typecode, k*m2->n+j);
                 sum += v1 * v2;
             }
-            outdata[j*m1->m+i] = sum;
+            outdata[i*m2->n+j] = sum;
         }
     }
     return MP_OBJ_FROM_PTR(out);
@@ -167,104 +169,8 @@ mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {
 
 MP_DEFINE_CONST_FUN_OBJ_2(linalg_dot_obj, linalg_dot);
 
-mp_obj_t linalg_zeros_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args, uint8_t kind) {
-    static const mp_arg_t allowed_args[] = {
-        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_obj = MP_OBJ_NULL} } ,
-        { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },
-    };
-
-    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);
-    
-    uint8_t dtype = args[1].u_int;
-    if(!MP_OBJ_IS_INT(args[0].u_obj) && !MP_OBJ_IS_TYPE(args[0].u_obj, &mp_type_tuple)) {
-        mp_raise_TypeError(translate("input argument must be an integer or a 2-tuple"));
-    }
-    ndarray_obj_t *ndarray = NULL;
-    if(MP_OBJ_IS_INT(args[0].u_obj)) {
-        size_t n = mp_obj_get_int(args[0].u_obj);
-        ndarray = create_new_ndarray(1, n, dtype);
-    } else if(MP_OBJ_IS_TYPE(args[0].u_obj, &mp_type_tuple)) {
-        mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(args[0].u_obj);
-        if(tuple->len != 2) {
-            mp_raise_TypeError(translate("input argument must be an integer or a 2-tuple"));
-        }
-        ndarray = create_new_ndarray(mp_obj_get_int(tuple->items[0]), 
-                                                  mp_obj_get_int(tuple->items[1]), dtype);
-    }
-    if(kind == 1) {
-        mp_obj_t one = mp_obj_new_int(1);
-        for(size_t i=0; i < ndarray->array->len; i++) {
-            mp_binary_set_val_array(dtype, ndarray->array->items, i, one);
-        }
-    }
-    return MP_OBJ_FROM_PTR(ndarray);
-}
-
-mp_obj_t linalg_zeros(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-    return linalg_zeros_ones(n_args, pos_args, kw_args, 0);
-}
-
-MP_DEFINE_CONST_FUN_OBJ_KW(linalg_zeros_obj, 0, linalg_zeros);
-
-mp_obj_t linalg_ones(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
-    return linalg_zeros_ones(n_args, pos_args, kw_args, 1);
-}
-
-MP_DEFINE_CONST_FUN_OBJ_KW(linalg_ones_obj, 0, linalg_ones);
-
-mp_obj_t linalg_eye(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_INT, {.u_int = 0} },
-        { MP_QSTR_M, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
-        { MP_QSTR_k, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0} },        
-        { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },
-    };
-
-    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);
-
-    size_t n = args[0].u_int, m;
-    int16_t k = args[2].u_int;
-    uint8_t dtype = args[3].u_int;
-    if(args[1].u_rom_obj == mp_const_none) {
-        m = n;
-    } else {
-        m = mp_obj_get_int(args[1].u_rom_obj);
-    }
-    
-    ndarray_obj_t *ndarray = create_new_ndarray(m, n, dtype);
-    mp_obj_t one = mp_obj_new_int(1);
-    size_t i = 0;
-    if((k >= 0) && (k < n)) {
-        while(k < n) {
-            mp_binary_set_val_array(dtype, ndarray->array->items, i*n+k, one);
-            k++;
-            i++;
-        }
-    } else if((k < 0) && (-k < m)) {
-        k = -k;
-        i = 0;
-        while(k < m) {
-            mp_binary_set_val_array(dtype, ndarray->array->items, k*n+i, one);
-            k++;
-            i++;
-        }
-    }
-    return MP_OBJ_FROM_PTR(ndarray);
-}
-
-MP_DEFINE_CONST_FUN_OBJ_KW(linalg_eye_obj, 0, linalg_eye);
-
-mp_obj_t linalg_det(mp_obj_t oin) {
-    if(!MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {
-        mp_raise_TypeError(translate("function defined for ndarrays only"));
-    }
-    ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
-    if(in->m != in->n) {
-        mp_raise_ValueError(translate("input must be square matrix"));
-    }
-    
+static mp_obj_t linalg_det(mp_obj_t oin) {
+    ndarray_obj_t *in = linalg_object_is_square(oin);  
     mp_float_t *tmp = m_new(mp_float_t, in->n*in->n);
     for(size_t i=0; i < in->array->len; i++){
         tmp[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);
@@ -295,14 +201,8 @@ mp_obj_t linalg_det(mp_obj_t oin) {
 
 MP_DEFINE_CONST_FUN_OBJ_1(linalg_det_obj, linalg_det);
 
-mp_obj_t linalg_eig(mp_obj_t oin) {
-    if(!MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {
-        mp_raise_TypeError(translate("function defined for ndarrays only"));
-    }
-    ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
-    if(in->m != in->n) {
-        mp_raise_ValueError(translate("input must be square matrix"));
-    }
+static mp_obj_t linalg_eig(mp_obj_t oin) {
+	ndarray_obj_t *in = linalg_object_is_square(oin);
     mp_float_t *array = m_new(mp_float_t, in->array->len);
     for(size_t i=0; i < in->array->len; i++) {
         array[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);
@@ -424,17 +324,80 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
 
 MP_DEFINE_CONST_FUN_OBJ_1(linalg_eig_obj, linalg_eig);
 
-#if !CIRCUITPY
+static mp_obj_t linalg_cholesky(mp_obj_t oin) {
+	ndarray_obj_t *in = MP_OBJ_TO_PTR(oin);
+	ndarray_obj_t *L = create_new_ndarray(in->n, in->n, NDARRAY_FLOAT);
+    mp_float_t *array = (mp_float_t *)L->array->items;
+
+	size_t pos = 0;
+    for(size_t m=0; m < in->m; m++) { // rows
+		for(size_t n=0; n < in->n; n++) { // columns
+			array[m*in->m+n] = ndarray_get_float_value(in->array->items, in->array->typecode, pos++);
+		}
+	}
+
+	// make sure the matrix is symmetric
+    for(size_t m=0; m < in->m; m++) { // rows
+        for(size_t n=m+1; n < in->n; n++) { // columns
+            // compare entry (m, n) to (n, m)
+            if(epsilon < MICROPY_FLOAT_C_FUN(fabs)(array[m*in->n + n] - array[n*in->n + m])) {
+                mp_raise_ValueError(translate("input matrix is asymmetric"));
+            }
+        }
+    }
+	
+	// this is actually not needed, but Cholesky in numpy returns the lower triangular matrix
+	for(size_t i=0; i < in->m; i++) { // rows
+		for(size_t j=i+1; j < in->n; j++) { // columns
+			array[i*in->m + j] = 0.0;
+		}
+	}	
+	mp_float_t sum = 0.0;
+	for(size_t i=0; i < in->m; i++) { // rows
+		for(size_t j=0; j <= i; j++) { // columns
+			sum = array[i*in->m + j];
+			for(size_t k=0; k < j; k++) {
+				sum -= array[i*in->n + k] * array[j*in->n + k];
+			}
+			if(i == j) {
+				if(sum <= 0.0) {
+					mp_raise_ValueError(translate("matrix is not positive definite"));
+				} else {
+					array[i*in->m+i] = MICROPY_FLOAT_C_FUN(sqrt)(sum);
+				}
+			} else {
+				array[i*in->m + j] = sum / array[j*in->m+j];
+			}
+		}
+	}
+	return MP_OBJ_FROM_PTR(L);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_1(linalg_cholesky_obj, linalg_cholesky);
+
+static mp_obj_t linalg_trace(mp_obj_t oin) {
+	ndarray_obj_t *ndarray = linalg_object_is_square(oin);
+	mp_float_t trace = 0.0;
+	for(size_t pos=0; pos < ndarray->array->len; pos+=(ndarray->m+1)) {
+		trace += ndarray_get_float_value(ndarray->array->items, ndarray->array->typecode, pos);
+	}
+	if(ndarray->array->typecode == NDARRAY_FLOAT) {
+		return mp_obj_new_float(trace);
+	}
+	return mp_obj_new_int_from_float(trace);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_1(linalg_trace_obj, linalg_trace);
+
 STATIC const mp_rom_map_elem_t ulab_linalg_globals_table[] = {
-    { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_linalg) },
-    { MP_ROM_QSTR(MP_QSTR_size), (mp_obj_t)&linalg_size_obj },
-    { MP_ROM_QSTR(MP_QSTR_inv), (mp_obj_t)&linalg_inv_obj },
-    { MP_ROM_QSTR(MP_QSTR_dot), (mp_obj_t)&linalg_dot_obj },
-    { MP_ROM_QSTR(MP_QSTR_zeros), (mp_obj_t)&linalg_zeros_obj },
-    { MP_ROM_QSTR(MP_QSTR_ones), (mp_obj_t)&linalg_ones_obj },
-    { MP_ROM_QSTR(MP_QSTR_eye), (mp_obj_t)&linalg_eye_obj },
-    { MP_ROM_QSTR(MP_QSTR_det), (mp_obj_t)&linalg_det_obj },
-    { MP_ROM_QSTR(MP_QSTR_eig), (mp_obj_t)&linalg_eig_obj },    
+	{ MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_linalg) },
+	{ MP_ROM_QSTR(MP_QSTR_size), (mp_obj_t)&linalg_size_obj },
+	{ MP_ROM_QSTR(MP_QSTR_inv), (mp_obj_t)&linalg_inv_obj },
+	{ MP_ROM_QSTR(MP_QSTR_dot), (mp_obj_t)&linalg_dot_obj },
+	{ MP_ROM_QSTR(MP_QSTR_det), (mp_obj_t)&linalg_det_obj },
+	{ MP_ROM_QSTR(MP_QSTR_eig), (mp_obj_t)&linalg_eig_obj },
+	{ MP_ROM_QSTR(MP_QSTR_cholesky), (mp_obj_t)&linalg_cholesky_obj },
+	{ MP_ROM_QSTR(MP_QSTR_trace), (mp_obj_t)&linalg_trace_obj },	
 };
 
 STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_linalg_globals, ulab_linalg_globals_table);
@@ -443,6 +406,5 @@ mp_obj_module_t ulab_linalg_module = {
     .base = { &mp_type_module },
     .globals = (mp_obj_dict_t*)&mp_module_ulab_linalg_globals,
 };
-#endif
 
 #endif

code/linalg.h

@@ -31,5 +31,11 @@ bool linalg_invert_matrix(mp_float_t *, size_t );
 
 extern mp_obj_module_t ulab_linalg_module;
 
+MP_DECLARE_CONST_FUN_OBJ_KW(linalg_size_obj);
+MP_DECLARE_CONST_FUN_OBJ_1(linalg_inv_obj);
+MP_DECLARE_CONST_FUN_OBJ_2(linalg_dot_obj);
+MP_DECLARE_CONST_FUN_OBJ_1(linalg_det_obj);
+MP_DECLARE_CONST_FUN_OBJ_1(linalg_eig_obj);
+
 #endif
 #endif

code/micropython.mk

@@ -3,12 +3,15 @@ USERMODULES_DIR := $(USERMOD_DIR)
 
 # Add all C files to SRC_USERMOD.
 SRC_USERMOD += $(USERMODULES_DIR)/ndarray.c
+SRC_USERMOD += $(USERMODULES_DIR)/create.c
 SRC_USERMOD += $(USERMODULES_DIR)/linalg.c
 SRC_USERMOD += $(USERMODULES_DIR)/vectorise.c
 SRC_USERMOD += $(USERMODULES_DIR)/poly.c
 SRC_USERMOD += $(USERMODULES_DIR)/fft.c
 SRC_USERMOD += $(USERMODULES_DIR)/numerical.c
 SRC_USERMOD += $(USERMODULES_DIR)/filter.c
+SRC_USERMOD += $(USERMODULES_DIR)/compare.c
+SRC_USERMOD += $(USERMODULES_DIR)/approx.c
 SRC_USERMOD += $(USERMODULES_DIR)/extras.c
 SRC_USERMOD += $(USERMODULES_DIR)/ulab.c
 

code/ndarray.c

@@ -9,6 +9,7 @@
  * Copyright (c) 2019-2020 Zoltán Vörös
 */
 
+#include <unistd.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -19,6 +20,9 @@
 #include "py/objtuple.h"
 #include "ndarray.h"
 
+mp_uint_t ndarray_print_threshold = NDARRAY_PRINT_THRESHOLD;
+mp_uint_t ndarray_print_edgeitems = NDARRAY_PRINT_EDGEITEMS;
+
 // This function is copied verbatim from objarray.c
 STATIC mp_obj_array_t *array_new(char typecode, size_t n) {
     int typecode_size = mp_binary_get_size('@', typecode, NULL);
@@ -52,6 +56,16 @@ mp_float_t ndarray_get_float_value(void *data, uint8_t typecode, size_t index) {
     }
 }
 
+bool ndarray_object_is_nditerable(mp_obj_t o_in) {
+    if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type) || 
+      MP_OBJ_IS_TYPE(o_in, &mp_type_tuple) || 
+      MP_OBJ_IS_TYPE(o_in, &mp_type_list) || 
+      MP_OBJ_IS_TYPE(o_in, &mp_type_range)) {
+        return true;
+    }
+    return false;
+}
+
 void fill_array_iterable(mp_float_t *array, mp_obj_t iterable) {
     mp_obj_iter_buf_t x_buf;
     mp_obj_t x_item, x_iterable = mp_getiter(iterable, &x_buf);
@@ -62,26 +76,53 @@ void fill_array_iterable(mp_float_t *array, mp_obj_t iterable) {
     }
 }
 
+mp_obj_t ndarray_set_printoptions(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+	static const mp_arg_t allowed_args[] = {
+		{ MP_QSTR_threshold, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} },
+		{ MP_QSTR_edgeitems, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} },
+	};
+
+	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);
+	if(args[0].u_rom_obj != mp_const_none) {
+		ndarray_print_threshold = mp_obj_get_int(args[0].u_rom_obj);
+	}
+	if(args[1].u_rom_obj != mp_const_none) {
+		ndarray_print_edgeitems = mp_obj_get_int(args[1].u_rom_obj);
+	}
+	return mp_const_none;
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(ndarray_set_printoptions_obj, 0, ndarray_set_printoptions);
+
+mp_obj_t ndarray_get_printoptions(void) {
+	mp_obj_t dict = mp_obj_new_dict(2);
+	mp_obj_dict_store(MP_OBJ_FROM_PTR(dict), MP_OBJ_NEW_QSTR(MP_QSTR_threshold), mp_obj_new_int(ndarray_print_threshold));
+	mp_obj_dict_store(MP_OBJ_FROM_PTR(dict), MP_OBJ_NEW_QSTR(MP_QSTR_edgeitems), mp_obj_new_int(ndarray_print_edgeitems));
+    return dict;
+}
+
+MP_DEFINE_CONST_FUN_OBJ_0(ndarray_get_printoptions_obj, ndarray_get_printoptions);
+
 void ndarray_print_row(const mp_print_t *print, mp_obj_array_t *data, size_t n0, size_t n) {
     mp_print_str(print, "[");
-    size_t i;
-    if(n < PRINT_MAX) { // if the array is short, print everything
+    if((n <= ndarray_print_threshold) || (n <= 2*ndarray_print_edgeitems)) { // if the array is short, print everything
         mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0), PRINT_REPR);
-        for(i=1; i<n; i++) {
+        for(size_t i=1; i < n; i++) {
             mp_print_str(print, ", ");
             mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+i), PRINT_REPR);
         }
     } else {
         mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0), PRINT_REPR);
-        for(i=1; i<3; i++) {
+        for(size_t i=1; i < ndarray_print_edgeitems; i++) {
             mp_print_str(print, ", ");
             mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+i), PRINT_REPR);
         }
         mp_printf(print, ", ..., ");
-        mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+n-3), PRINT_REPR);
-        for(size_t i=1; i<3; i++) {
+        mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+n-ndarray_print_edgeitems), PRINT_REPR);
+        for(size_t i=1; i < ndarray_print_edgeitems; i++) {
             mp_print_str(print, ", ");
-            mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+n-3+i), PRINT_REPR);
+            mp_obj_print_helper(print, mp_binary_get_val_array(data->typecode, data->items, n0+n-ndarray_print_edgeitems+i), PRINT_REPR);
         }
     }
     mp_print_str(print, "]");
@@ -98,7 +139,7 @@ void ndarray_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t ki
         if((self->m == 1) || (self->n == 1)) {
             ndarray_print_row(print, self->array, 0, self->array->len);
         } else {
-            // TODO: add vertical ellipses for the case, when self->m > PRINT_MAX
+            // TODO: add vertical ellipses
             mp_print_str(print, "[");
             ndarray_print_row(print, self->array, 0, self->n);
             for(size_t i=1; i < self->m; i++) {
@@ -159,16 +200,38 @@ STATIC uint8_t ndarray_init_helper(size_t n_args, const mp_obj_t *pos_args, mp_m
     };
     
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
     
     uint8_t dtype = args[1].u_int;
     return dtype;
 }
 
-STATIC mp_obj_t ndarray_make_new_core(const mp_obj_type_t *type, size_t n_args, size_t n_kw, const mp_obj_t *args, mp_map_t *kw_args) {
+STATIC mp_obj_t ndarray_make_new_core(size_t n_args, const mp_obj_t *args, mp_map_t *kw_args) {
     uint8_t dtype = ndarray_init_helper(n_args, args, kw_args);
+    
+	if(MP_OBJ_IS_TYPE(args[0], &ulab_ndarray_type)) {
+		ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0]);
+		if(dtype == ndarray->array->typecode) {
+			return ndarray_copy(args[0]);
+		}
+		ndarray_obj_t *ndarray_new = create_new_ndarray(ndarray->m, ndarray->n, dtype);
+		mp_obj_t item;
+		if((ndarray->array->typecode == NDARRAY_FLOAT) &&(dtype != NDARRAY_FLOAT)) {
+			for(size_t i=0; i < ndarray->array->len; i++) {
+				mp_float_t f = ndarray_get_float_value(ndarray->array->items, ndarray->array->typecode, i);
+				item = mp_obj_new_int((int32_t)MICROPY_FLOAT_C_FUN(floor)(f));
+				mp_binary_set_val_array(dtype, ndarray_new->array->items, i, item);
+			}
+		} else {
+			for(size_t i=0; i < ndarray->array->len; i++) {
+				item = mp_binary_get_val_array(ndarray->array->typecode, ndarray->array->items, i);
+				mp_binary_set_val_array(dtype, ndarray_new->array->items, i, item);
+			}
+		}
+		return MP_OBJ_FROM_PTR(ndarray_new);
+	}
 
-    size_t len1, len2=0, i=0;
+    size_t len1, len2 = 0, i = 0;
     mp_obj_t len_in = mp_obj_len_maybe(args[0]);
     if (len_in == MP_OBJ_NULL) {
         mp_raise_ValueError(translate("first argument must be an iterable"));
@@ -186,7 +249,7 @@ STATIC mp_obj_t ndarray_make_new_core(const mp_obj_type_t *type, size_t n_args,
         if(len_in != MP_OBJ_NULL) { // indeed, this seems to be an iterable
             // Next, we have to check, whether all elements in the outer loop have the same length
             if(i > 0) {
-                if(len2 != MP_OBJ_SMALL_INT_VALUE(len_in)) {
+                if(len2 != (size_t)MP_OBJ_SMALL_INT_VALUE(len_in)) {
                     mp_raise_ValueError(translate("iterables are not of the same length"));
                 }
             }
@@ -214,42 +277,44 @@ STATIC mp_obj_t ndarray_make_new_core(const mp_obj_type_t *type, size_t n_args,
 
 #ifdef CIRCUITPY
 mp_obj_t ndarray_make_new(const mp_obj_type_t *type, size_t n_args, const mp_obj_t *args, mp_map_t *kw_args) {
+    (void) type;
     mp_arg_check_num(n_args, kw_args, 1, 2, true);
     size_t n_kw = 0;
     if (kw_args != 0) {
         n_kw = kw_args->used;
     }
     mp_map_init_fixed_table(kw_args, n_kw, args + n_args);
-    return ndarray_make_new_core(type, n_args, n_kw, args, kw_args);
+    return ndarray_make_new_core(n_args, args, kw_args);
 }
 #else
 mp_obj_t ndarray_make_new(const mp_obj_type_t *type, size_t n_args, size_t n_kw, const mp_obj_t *args) {
+    (void) type;
     mp_arg_check_num(n_args, n_kw, 1, 2, true);
     mp_map_t kw_args;
     mp_map_init_fixed_table(&kw_args, n_kw, args + n_args);
-    return ndarray_make_new_core(type, n_args, n_kw, args, &kw_args);
+    return ndarray_make_new_core(n_args, args, &kw_args);
 }
 #endif
 
-size_t slice_length(mp_bound_slice_t slice) {
-    int32_t len, correction = 1;
+static size_t slice_length(mp_bound_slice_t slice) {
+    ssize_t len, correction = 1;
     if(slice.step > 0) correction = -1;
-    len = (slice.stop - slice.start + (slice.step + correction)) / slice.step;
+    len = (ssize_t)(slice.stop - slice.start + (slice.step + correction)) / slice.step;
     if(len < 0) return 0;
     return (size_t)len;
 }
 
-size_t true_length(mp_obj_t bool_list) {
+static size_t true_length(mp_obj_t bool_list) {
     // returns the number of Trues in a Boolean list
     // I wonder, wouldn't this be faster, if we looped through bool_list->items instead?
     mp_obj_iter_buf_t iter_buf;
     mp_obj_t item, iterable = mp_getiter(bool_list, &iter_buf);
     size_t trues = 0;
     while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
-        if(!MP_OBJ_IS_TYPE(item, &mp_type_bool)) {
-            // numpy seems to be a little bit inconsistent in when an index is considered
-            // to be True/False. Bail out immediately, if the items are not True/False
-            return 0;
+        if(!mp_obj_is_bool(item)) {
+                // numpy seems to be a little bit inconsistent in when an index is considered
+                // to be True/False. Bail out immediately, if the items are not True/False
+                mp_raise_TypeError(translate("wrong index type"));
         }
         if(mp_obj_is_true(item)) {
             trues++;
@@ -258,13 +323,13 @@ size_t true_length(mp_obj_t bool_list) {
     return trues;
 }
 
-mp_bound_slice_t generate_slice(mp_uint_t n, mp_obj_t index) {
+static mp_bound_slice_t generate_slice(mp_int_t n, mp_obj_t index) {
     // micropython seems to have difficulties with negative steps
     mp_bound_slice_t slice;
     if(MP_OBJ_IS_TYPE(index, &mp_type_slice)) {
-        mp_seq_get_fast_slice_indexes(n, index, &slice);
+        mp_obj_slice_indices(index, n, &slice);
     } else if(MP_OBJ_IS_INT(index)) {
-        int32_t _index = mp_obj_get_int(index);
+        mp_int_t _index = mp_obj_get_int(index);
         if(_index < 0) {
             _index += n;
         } 
@@ -280,7 +345,7 @@ mp_bound_slice_t generate_slice(mp_uint_t n, mp_obj_t index) {
     return slice;
 }
 
-mp_bound_slice_t simple_slice(int16_t start, int16_t stop, int16_t step) {
+static mp_bound_slice_t simple_slice(int16_t start, int16_t stop, int16_t step) {
     mp_bound_slice_t slice;
     slice.start = start;
     slice.stop = stop;
@@ -288,18 +353,18 @@ mp_bound_slice_t simple_slice(int16_t start, int16_t stop, int16_t step) {
     return slice;
 }
 
-void insert_binary_value(ndarray_obj_t *ndarray, size_t nd_index, ndarray_obj_t *values, size_t value_index) {
+static void insert_binary_value(ndarray_obj_t *ndarray, size_t nd_index, ndarray_obj_t *values, size_t value_index) {
     // there is probably a more elegant implementation...
     mp_obj_t tmp = mp_binary_get_val_array(values->array->typecode, values->array->items, value_index);
     if((values->array->typecode == NDARRAY_FLOAT) && (ndarray->array->typecode != NDARRAY_FLOAT)) {
         // workaround: rounding seems not to work in the arm compiler
-        int32_t x = (int32_t)floorf(mp_obj_get_float(tmp)+0.5);
+        int32_t x = (int32_t)MICROPY_FLOAT_C_FUN(floor)(mp_obj_get_float(tmp)+0.5);
         tmp = mp_obj_new_int(x);
     }
     mp_binary_set_val_array(ndarray->array->typecode, ndarray->array->items, nd_index, tmp); 
 }
 
-mp_obj_t insert_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n, 
+static mp_obj_t insert_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n, 
                             mp_bound_slice_t row, mp_bound_slice_t column, 
                             mp_obj_t row_list, mp_obj_t column_list, 
                             ndarray_obj_t *values) {
@@ -367,7 +432,8 @@ mp_obj_t insert_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n,
         } else { // columns are indexed by a list
             mp_obj_iter_buf_t column_iter_buf;
             mp_obj_t column_item, column_iterable;
-            size_t j = 0, cindex = 0;
+            size_t j = 0;
+            cindex = 0;
             while((row_item = mp_iternext(row_iterable)) != MP_OBJ_STOP_ITERATION) {
                 if(mp_obj_is_true(row_item)) {
                     column_iterable = mp_getiter(column_list, &column_iter_buf);                   
@@ -387,14 +453,10 @@ mp_obj_t insert_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n,
     return mp_const_none;
 }
 
-mp_obj_t iterate_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n, 
+static mp_obj_t iterate_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n, 
                             mp_bound_slice_t row, mp_bound_slice_t column, 
                             mp_obj_t row_list, mp_obj_t column_list, 
                             ndarray_obj_t *values) {
-    if((m == 0) || (n == 0)) {
-        mp_raise_msg(&mp_type_IndexError, translate("empty index range"));
-    }
-
     if(values != NULL) {
         return insert_slice_list(ndarray, m, n, row, column, row_list, column_list, values);
     }
@@ -453,7 +515,8 @@ mp_obj_t iterate_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n,
         } else { // columns are indexed by a list
             mp_obj_iter_buf_t column_iter_buf;
             mp_obj_t column_item, column_iterable;
-            size_t j = 0, cindex = 0;
+            size_t j = 0;
+            cindex = 0;
             while((row_item = mp_iternext(row_iterable)) != MP_OBJ_STOP_ITERATION) {
                 if(mp_obj_is_true(row_item)) {
                     column_iterable = mp_getiter(column_list, &column_iter_buf);                   
@@ -473,7 +536,7 @@ mp_obj_t iterate_slice_list(ndarray_obj_t *ndarray, size_t m, size_t n,
     return MP_OBJ_FROM_PTR(out);
 }
 
-mp_obj_t ndarray_get_slice(ndarray_obj_t *ndarray, mp_obj_t index, ndarray_obj_t *values) {
+static mp_obj_t ndarray_get_slice(ndarray_obj_t *ndarray, mp_obj_t index, ndarray_obj_t *values) {
     mp_bound_slice_t row_slice = simple_slice(0, 0, 1), column_slice = simple_slice(0, 0, 1);
 
     size_t m = 0, n = 0;
@@ -590,7 +653,7 @@ typedef struct _mp_obj_ndarray_it_t {
     size_t cur;
 } mp_obj_ndarray_it_t;
 
-mp_obj_t ndarray_iternext(mp_obj_t self_in) {
+static mp_obj_t ndarray_iternext(mp_obj_t self_in) {
     mp_obj_ndarray_it_t *self = MP_OBJ_TO_PTR(self_in);
     ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(self->ndarray);
     // TODO: in numpy, ndarrays are iterated with respect to the first axis. 
@@ -686,155 +749,176 @@ mp_obj_t ndarray_flatten(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a
 }
 
 // Binary operations
+ndarray_obj_t *ndarray_from_mp_obj(mp_obj_t obj) {
+	// creates an ndarray from an micropython int or float
+	// if the input is an ndarray, it is returned
+	ndarray_obj_t *ndarray;
+	if(MP_OBJ_IS_INT(obj)) {
+		int32_t ivalue = mp_obj_get_int(obj);
+		if((ivalue > 0) && (ivalue < 256)) {
+			CREATE_SINGLE_ITEM(ndarray, uint8_t, NDARRAY_UINT8, ivalue);
+		} else if((ivalue > 255) && (ivalue < 65535)) {
+			CREATE_SINGLE_ITEM(ndarray, uint16_t, NDARRAY_UINT16, ivalue);
+		} else if((ivalue < 0) && (ivalue > -128)) {
+			CREATE_SINGLE_ITEM(ndarray, int8_t, NDARRAY_INT8, ivalue);
+		} else if((ivalue < -127) && (ivalue > -32767)) {
+			CREATE_SINGLE_ITEM(ndarray, int16_t, NDARRAY_INT16, ivalue);
+		} else { // the integer value clearly does not fit the ulab types, so move on to float
+			CREATE_SINGLE_ITEM(ndarray, mp_float_t, NDARRAY_FLOAT, ivalue);
+		}
+	} else if(mp_obj_is_float(obj)) {
+		mp_float_t fvalue = mp_obj_get_float(obj);
+		CREATE_SINGLE_ITEM(ndarray, mp_float_t, NDARRAY_FLOAT, fvalue);
+	} else if(MP_OBJ_IS_TYPE(obj, &ulab_ndarray_type)){
+		ndarray = MP_OBJ_TO_PTR(obj);
+	} else {
+		mp_raise_TypeError(translate("wrong operand type"));
+	}
+	return ndarray;
+}
 
-mp_obj_t ndarray_binary_op(mp_binary_op_t op, mp_obj_t lhs, mp_obj_t rhs) {
-//    if(op == MP_BINARY_OP_REVERSE_ADD) {
- //       return ndarray_binary_op(MP_BINARY_OP_ADD, rhs, lhs);
-  //  }    
+mp_obj_t ndarray_binary_op(mp_binary_op_t _op, mp_obj_t lhs, mp_obj_t rhs) {
+	// if the ndarray stands on the right hand side of the expression, simply swap the operands
+	ndarray_obj_t *ol, *or;
+	mp_binary_op_t op = _op;
+	if((op == MP_BINARY_OP_REVERSE_ADD) || (op == MP_BINARY_OP_REVERSE_MULTIPLY) || 
+		(op == MP_BINARY_OP_REVERSE_POWER) || (op == MP_BINARY_OP_REVERSE_SUBTRACT) || 
+		(op == MP_BINARY_OP_REVERSE_TRUE_DIVIDE)) {
+		ol = ndarray_from_mp_obj(rhs);
+		or = ndarray_from_mp_obj(lhs);		
+	} else {
+		ol = ndarray_from_mp_obj(lhs);
+		or = ndarray_from_mp_obj(rhs);
+	}
+	if(op == MP_BINARY_OP_REVERSE_ADD) {
+		op = MP_BINARY_OP_ADD;
+	} else if(op == MP_BINARY_OP_REVERSE_MULTIPLY) {
+		op = MP_BINARY_OP_MULTIPLY;
+	} else if(op == MP_BINARY_OP_REVERSE_POWER) {
+		op = MP_BINARY_OP_POWER;
+	} else if(op == MP_BINARY_OP_REVERSE_SUBTRACT) {
+		op = MP_BINARY_OP_SUBTRACT;
+	} else if(op == MP_BINARY_OP_REVERSE_TRUE_DIVIDE) {
+		op = MP_BINARY_OP_TRUE_DIVIDE;
+	}
     // One of the operands is a scalar
     // TODO: conform to numpy with the upcasting
     // TODO: implement in-place operators
-    mp_obj_t RHS = MP_OBJ_NULL;
-    bool rhs_is_scalar = true;
-    if(MP_OBJ_IS_INT(rhs)) {
-        int32_t ivalue = mp_obj_get_int(rhs);
-        if((ivalue > 0) && (ivalue < 256)) {
-            CREATE_SINGLE_ITEM(RHS, uint8_t, NDARRAY_UINT8, ivalue);
-        } else if((ivalue > 255) && (ivalue < 65535)) {
-            CREATE_SINGLE_ITEM(RHS, uint16_t, NDARRAY_UINT16, ivalue);
-        } else if((ivalue < 0) && (ivalue > -128)) {
-            CREATE_SINGLE_ITEM(RHS, int8_t, NDARRAY_INT8, ivalue);
-        } else if((ivalue < -127) && (ivalue > -32767)) {
-            CREATE_SINGLE_ITEM(RHS, int16_t, NDARRAY_INT16, ivalue);
-        } else { // the integer value clearly does not fit the ulab types, so move on to float
-            CREATE_SINGLE_ITEM(RHS, mp_float_t, NDARRAY_FLOAT, ivalue);
-        }
-    } else if(mp_obj_is_float(rhs)) {
-        mp_float_t fvalue = mp_obj_get_float(rhs);        
-        CREATE_SINGLE_ITEM(RHS, mp_float_t, NDARRAY_FLOAT, fvalue);
-    } else {
-        RHS = rhs;
-        rhs_is_scalar = false;
-    }
-    //else 
-    if(MP_OBJ_IS_TYPE(lhs, &ulab_ndarray_type) && MP_OBJ_IS_TYPE(RHS, &ulab_ndarray_type)) { 
-        // next, the ndarray stuff
-        ndarray_obj_t *ol = MP_OBJ_TO_PTR(lhs);
-        ndarray_obj_t *or = MP_OBJ_TO_PTR(RHS);
-        if(!rhs_is_scalar && ((ol->m != or->m) || (ol->n != or->n))) {
+	// these are partial broadcasting rules: either the two arrays 
+	// are of the same shape, or one of them is of length 1
+	if(((ol->m != or->m) || (ol->n != or->n))) {
+		if((ol->array->len != 1) && (or->array->len != 1)) {
+			if(op == MP_BINARY_OP_EQUAL) {
+				return mp_const_false;
+			} else if(op == MP_BINARY_OP_NOT_EQUAL) {
+				return mp_const_true;
+			}
             mp_raise_ValueError(translate("operands could not be broadcast together"));
-        }
-        // At this point, the operands should have the same shape
-        switch(op) {
-            case MP_BINARY_OP_EQUAL:
-                // Two arrays are equal, if their shape, typecode, and elements are equal
-                if((ol->m != or->m) || (ol->n != or->n) || (ol->array->typecode != or->array->typecode)) {
-                    return mp_const_false;
-                } else {
-                    size_t i = ol->bytes;
-                    uint8_t *l = (uint8_t *)ol->array->items;
-                    uint8_t *r = (uint8_t *)or->array->items;
-                    while(i) { // At this point, we can simply compare the bytes, the type is irrelevant
-                        if(*l++ != *r++) {
-                            return mp_const_false;
-                        }
-                        i--;
-                    }
-                    return mp_const_true;
-                }
-                break;
-            case MP_BINARY_OP_LESS:
-            case MP_BINARY_OP_LESS_EQUAL:
-            case MP_BINARY_OP_MORE:
-            case MP_BINARY_OP_MORE_EQUAL:
-            case MP_BINARY_OP_ADD:
-            case MP_BINARY_OP_SUBTRACT:
-            case MP_BINARY_OP_TRUE_DIVIDE:
-            case MP_BINARY_OP_MULTIPLY:
-                // TODO: I believe, this part can be made significantly smaller (compiled size)
-                // by doing only the typecasting in the large ifs, and moving the loops outside
-                // These are the upcasting rules
-                // float always becomes float
-                // operation on identical types preserves type
-                // uint8 + int8 => int16
-                // uint8 + int16 => int16
-                // uint8 + uint16 => uint16
-                // int8 + int16 => int16
-                // int8 + uint16 => uint16
-                // uint16 + int16 => float
-                // The parameters of RUN_BINARY_LOOP are 
-                // typecode of result, type_out, type_left, type_right, lhs operand, rhs operand, operator
-                if(ol->array->typecode == NDARRAY_UINT8) {
-                    if(or->array->typecode == NDARRAY_UINT8) {
-                        RUN_BINARY_LOOP(NDARRAY_UINT8, uint8_t, uint8_t, uint8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT8) {
-                        RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_UINT16) {
-                        RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint8_t, uint16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT16) {
-                        RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_FLOAT) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op);
-                    }
-                } else if(ol->array->typecode == NDARRAY_INT8) {
-                    if(or->array->typecode == NDARRAY_UINT8) {
-                        RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT8) {
-                        RUN_BINARY_LOOP(NDARRAY_INT8, int8_t, int8_t, int8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_UINT16) {
-                        RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT16) {
-                        RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, int16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_FLOAT) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int8_t, mp_float_t, ol, or, op);
-                    }                
-                } else if(ol->array->typecode == NDARRAY_UINT16) {
-                    if(or->array->typecode == NDARRAY_UINT8) {
-                        RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT8) {
-                        RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, int8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_UINT16) {
-                        RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT16) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, int16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_FLOAT) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op);
-                    }
-                } else if(ol->array->typecode == NDARRAY_INT16) {
-                    if(or->array->typecode == NDARRAY_UINT8) {
-                        RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, uint8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT8) {
-                        RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_UINT16) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int16_t, uint16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT16) {
-                        RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_FLOAT) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, mp_float_t, ol, or, op);
-                    }
-                } else if(ol->array->typecode == NDARRAY_FLOAT) {
-                    if(or->array->typecode == NDARRAY_UINT8) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT8) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int8_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_UINT16) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_INT16) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int16_t, ol, or, op);
-                    } else if(or->array->typecode == NDARRAY_FLOAT) {
-                        RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, mp_float_t, ol, or, op);
-                    }
-                } else { // this should never happen
-                    mp_raise_TypeError(translate("wrong input type"));
-                }
-                // this instruction should never be reached, but we have to make the compiler happy
-                return MP_OBJ_NULL; 
-            default:
-                return MP_OBJ_NULL; // op not supported                                                        
-        }
-    } else {
-        mp_raise_TypeError(translate("wrong operand type on the right hand side"));
-    }
+		}
+	}
+	uint8_t linc = ol->array->len == 1 ? 0 : 1;
+	uint8_t rinc = or->array->len == 1 ? 0 : 1;
+	// do the partial broadcasting here
+	size_t m = MAX(ol->m, or->m);
+	size_t n = MAX(ol->n, or->n);
+	size_t len = MAX(ol->array->len, or->array->len);
+	if((ol->array->len == 0) || (or->array->len == 0)) {
+		len = 0;
+	}
+	switch(op) {
+		case MP_BINARY_OP_EQUAL:
+		case MP_BINARY_OP_NOT_EQUAL:
+		case MP_BINARY_OP_LESS:
+		case MP_BINARY_OP_LESS_EQUAL:
+		case MP_BINARY_OP_MORE:
+		case MP_BINARY_OP_MORE_EQUAL:
+		case MP_BINARY_OP_ADD:
+		case MP_BINARY_OP_SUBTRACT:
+		case MP_BINARY_OP_TRUE_DIVIDE:
+		case MP_BINARY_OP_MULTIPLY:
+		case MP_BINARY_OP_POWER:       
+			// TODO: I believe, this part can be made significantly smaller (compiled size)
+			// by doing only the typecasting in the large ifs, and moving the loops outside
+			// These are the upcasting rules
+			// float always becomes float
+			// operation on identical types preserves type
+			// uint8 + int8 => int16
+			// uint8 + int16 => int16
+			// uint8 + uint16 => uint16
+			// int8 + int16 => int16
+			// int8 + uint16 => uint16
+			// uint16 + int16 => float
+			// The parameters of RUN_BINARY_LOOP are 
+			// typecode of result, type_out, type_left, type_right, lhs operand, rhs operand, operator
+			if(ol->array->typecode == NDARRAY_UINT8) {
+				if(or->array->typecode == NDARRAY_UINT8) {
+					RUN_BINARY_LOOP(NDARRAY_UINT8, uint8_t, uint8_t, uint8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT8) {
+					RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_UINT16) {
+					RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint8_t, uint16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT16) {
+					RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, uint8_t, int16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_FLOAT) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint8_t, mp_float_t, ol, or, op, m, n, len, linc, rinc);
+				}
+			} else if(ol->array->typecode == NDARRAY_INT8) {
+				if(or->array->typecode == NDARRAY_UINT8) {
+					RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT8) {
+					RUN_BINARY_LOOP(NDARRAY_INT8, int8_t, int8_t, int8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_UINT16) {
+					RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, uint16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT16) {
+					RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int8_t, int16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_FLOAT) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int8_t, mp_float_t, ol, or, op, m, n, len, linc, rinc);
+				}                
+			} else if(ol->array->typecode == NDARRAY_UINT16) {
+				if(or->array->typecode == NDARRAY_UINT8) {
+					RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT8) {
+					RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, int8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_UINT16) {
+					RUN_BINARY_LOOP(NDARRAY_UINT16, uint16_t, uint16_t, uint16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT16) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, int16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_FLOAT) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, mp_float_t, ol, or, op, m, n, len, linc, rinc);
+				}
+			} else if(ol->array->typecode == NDARRAY_INT16) {
+				if(or->array->typecode == NDARRAY_UINT8) {
+					RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, uint8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT8) {
+					RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_UINT16) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, int16_t, uint16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT16) {
+					RUN_BINARY_LOOP(NDARRAY_INT16, int16_t, int16_t, int16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_FLOAT) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, uint16_t, mp_float_t, ol, or, op, m, n, len, linc, rinc);
+				}
+			} else if(ol->array->typecode == NDARRAY_FLOAT) {
+				if(or->array->typecode == NDARRAY_UINT8) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT8) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int8_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_UINT16) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, uint16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_INT16) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, int16_t, ol, or, op, m, n, len, linc, rinc);
+				} else if(or->array->typecode == NDARRAY_FLOAT) {
+					RUN_BINARY_LOOP(NDARRAY_FLOAT, mp_float_t, mp_float_t, mp_float_t, ol, or, op, m, n, len, linc, rinc);
+				}
+			}
+			// this instruction should never be reached, but we have to make the compiler happy
+			return MP_OBJ_NULL;
+			break;
+		default:
+			return MP_OBJ_NULL; // op not supported
+			break;
+	}
+	return MP_OBJ_NULL;
 }
 
 mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) {
@@ -856,9 +940,11 @@ mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) {
             // we can invert the content byte by byte, there is no need to distinguish 
             // between different typecodes
             ndarray = MP_OBJ_TO_PTR(ndarray_copy(self_in));
-            uint8_t *array = (uint8_t *)ndarray->array->items;
-            for(size_t i=0; i < self->bytes; i++) array[i] = ~array[i];
-            return MP_OBJ_FROM_PTR(ndarray);
+            {
+                uint8_t *array = (uint8_t *)ndarray->array->items;
+                for(size_t i=0; i < self->bytes; i++) array[i] = ~array[i];
+                return MP_OBJ_FROM_PTR(ndarray);
+            }
             break;
         
         case MP_UNARY_OP_NEGATIVE:
@@ -914,34 +1000,28 @@ mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) {
 
 mp_obj_t ndarray_transpose(mp_obj_t self_in) {
     ndarray_obj_t *self = MP_OBJ_TO_PTR(self_in);
-    // the size of a single item in the array
-    uint8_t _sizeof = mp_binary_get_size('@', self->array->typecode, NULL);
     
-    // NOTE: 
-    // if the matrices are square, we can simply swap items, but 
-    // generic matrices can't be transposed in place, so we have to 
-    // declare a temporary variable
-    
-    // NOTE: 
     //  In the old matrix, the coordinate (m, n) is m*self->n + n
     //  We have to assign this to the coordinate (n, m) in the new 
     //  matrix, i.e., to n*self->m + m (since the new matrix has self->m columns)
-    
+
+	ndarray_obj_t *ndarray = create_new_ndarray(self->n,  self->m, self->array->typecode);
     // one-dimensional arrays can be transposed by simply swapping the dimensions
-    if((self->m != 1) && (self->n != 1)) {
-        uint8_t *c = (uint8_t *)self->array->items;
-        // self->bytes is the size of the bytearray, irrespective of the typecode
-        uint8_t *tmp = m_new(uint8_t, self->bytes);
+    if((self->m == 1) || (self->n == 1)) {
+        memcpy(ndarray->array->items, self->array->items, self->bytes);
+	} else {
+		// the size of a single item in the array
+		uint8_t itemsize = mp_binary_get_size('@', self->array->typecode, NULL);
+        uint8_t *sarray = (uint8_t *)self->array->items;
+        uint8_t *narray = (uint8_t *)ndarray->array->items;
         for(size_t m=0; m < self->m; m++) {
             for(size_t n=0; n < self->n; n++) {
-                memcpy(tmp+_sizeof*(n*self->m + m), c+_sizeof*(m*self->n + n), _sizeof);
+                memcpy(narray+itemsize*(n*self->m + m), sarray, itemsize);
+                sarray += itemsize;
             }
         }
-        memcpy(self->array->items, tmp, self->bytes);
-        m_del(uint8_t, tmp, self->bytes);
     } 
-    SWAP(size_t, self->m, self->n);
-    return mp_const_none;
+    return MP_OBJ_FROM_PTR(ndarray);
 }
 
 MP_DEFINE_CONST_FUN_OBJ_1(ndarray_transpose_obj, ndarray_transpose);

code/ndarray.h

@@ -17,7 +17,8 @@
 #include "py/objstr.h"
 #include "py/objlist.h"
 
-#define PRINT_MAX  10
+#define NDARRAY_PRINT_THRESHOLD  	10
+#define NDARRAY_PRINT_EDGEITEMS		3
 
 #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
 #define FLOAT_TYPECODE 'f'
@@ -25,8 +26,11 @@
 #define FLOAT_TYPECODE 'd'
 #endif
 
-#if !CIRCUITPY
-#define translate(x) x
+#if CIRCUITPY
+#define mp_obj_is_bool(o) (MP_OBJ_IS_TYPE((o), &mp_type_bool))
+#define mp_obj_is_int(x) (MP_OBJ_IS_INT((x)))
+#else
+#define translate(x) MP_ERROR_TEXT(x)
 #endif
 
 #define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }
@@ -52,8 +56,13 @@ typedef struct _ndarray_obj_t {
 mp_obj_t mp_obj_new_ndarray_iterator(mp_obj_t , size_t , mp_obj_iter_buf_t *);
 
 mp_float_t ndarray_get_float_value(void *, uint8_t , size_t );
+bool ndarray_object_is_nditerable(mp_obj_t );
 void fill_array_iterable(mp_float_t *, mp_obj_t );
 
+mp_obj_t ndarray_set_printoptions(size_t , const mp_obj_t *, mp_map_t *);
+MP_DECLARE_CONST_FUN_OBJ_KW(ndarray_set_printoptions_obj);
+mp_obj_t ndarray_get_printoptions(void);
+MP_DECLARE_CONST_FUN_OBJ_0(ndarray_get_printoptions_obj);
 void ndarray_print_row(const mp_print_t *, mp_obj_array_t *, size_t , size_t );
 void ndarray_print(const mp_print_t *, mp_obj_t , mp_print_kind_t );
 void ndarray_assign_elements(mp_obj_array_t *, mp_obj_t , uint8_t , size_t *);
@@ -84,6 +93,7 @@ MP_DECLARE_CONST_FUN_OBJ_1(ndarray_transpose_obj);
 mp_int_t ndarray_get_buffer(mp_obj_t obj, mp_buffer_info_t *bufinfo, mp_uint_t flags);
 //void ndarray_attributes(mp_obj_t , qstr , mp_obj_t *);
 
+ndarray_obj_t *ndarray_from_mp_obj(mp_obj_t );
 
 #define CREATE_SINGLE_ITEM(outarray, type, typecode, value) do {\
     ndarray_obj_t *tmp = create_new_ndarray(1, 1, (typecode));\
@@ -99,45 +109,38 @@ mp_int_t ndarray_get_buffer(mp_obj_t obj, mp_buffer_info_t *bufinfo, mp_uint_t f
     should work outside the loop, but it doesn't. Go figure! 
 */
 
-#define RUN_BINARY_LOOP(typecode, type_out, type_left, type_right, ol, or, op) do {\
+#define RUN_BINARY_LOOP(typecode, type_out, type_left, type_right, ol, or, op, m, n, len, linc, rinc) do {\
     type_left *left = (type_left *)(ol)->array->items;\
     type_right *right = (type_right *)(or)->array->items;\
-    uint8_t inc = 0;\
-    if((or)->array->len > 1) inc = 1;\
-    if(((op) == MP_BINARY_OP_ADD) || ((op) == MP_BINARY_OP_SUBTRACT) || ((op) == MP_BINARY_OP_MULTIPLY)) {\
-        ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, typecode);\
+    if(((op) == MP_BINARY_OP_ADD) || ((op) == MP_BINARY_OP_SUBTRACT) || ((op) == MP_BINARY_OP_MULTIPLY) || ((op) == MP_BINARY_OP_POWER)) {\
+        ndarray_obj_t *out = create_new_ndarray((m), (n), (typecode));\
         type_out *(odata) = (type_out *)out->array->items;\
-        if((op) == MP_BINARY_OP_ADD) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] + right[j];}\
-        if((op) == MP_BINARY_OP_SUBTRACT) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] - right[j];}\
-        if((op) == MP_BINARY_OP_MULTIPLY) { for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = left[i] * right[j];}\
+        if((op) == MP_BINARY_OP_ADD) { for(size_t i=0; i < (len); i++, left+=linc, right+=rinc) *odata++ = *left + *right; }\
+		else if((op) == MP_BINARY_OP_MULTIPLY) { for(size_t i=0; i < (len); i++, left+=linc, right+=rinc) *odata++ = *left * *right; }\
+        else if((op) == MP_BINARY_OP_POWER) { for(size_t i=0; i < (len); i++, left+=linc, right+=rinc) *odata++ = (type_out)MICROPY_FLOAT_C_FUN(pow)(*left, *right); }\
+        else if((op) == MP_BINARY_OP_SUBTRACT) { for(size_t i=0; i < (len); i++, left+=linc, right+=rinc) *odata++ = *left - *right; }\
         return MP_OBJ_FROM_PTR(out);\
-    } else if((op) == MP_BINARY_OP_TRUE_DIVIDE) {\
-        ndarray_obj_t *out = create_new_ndarray(ol->m, ol->n, NDARRAY_FLOAT);\
+	} else if((op) == MP_BINARY_OP_TRUE_DIVIDE) {\
+        ndarray_obj_t *out = create_new_ndarray((m), (n), NDARRAY_FLOAT);\
         mp_float_t *odata = (mp_float_t *)out->array->items;\
-        for(size_t i=0, j=0; i < (ol)->array->len; i++, j+=inc) odata[i] = (mp_float_t)left[i]/(mp_float_t)right[j];\
+        for(size_t i=0; i < (len); i++, left+=linc, right+=rinc) {*odata++ = (mp_float_t)(*left)/(mp_float_t)(*right);}\
         return MP_OBJ_FROM_PTR(out);\
-    } else if(((op) == MP_BINARY_OP_LESS) || ((op) == MP_BINARY_OP_LESS_EQUAL) ||  \
-             ((op) == MP_BINARY_OP_MORE) || ((op) == MP_BINARY_OP_MORE_EQUAL)) {\
+	} else if(((op) == MP_BINARY_OP_LESS) || ((op) == MP_BINARY_OP_LESS_EQUAL) ||  \
+             ((op) == MP_BINARY_OP_MORE) || ((op) == MP_BINARY_OP_MORE_EQUAL) || \
+             ((op) == MP_BINARY_OP_EQUAL) || ((op) == MP_BINARY_OP_NOT_EQUAL)) {\
         mp_obj_t out_list = mp_obj_new_list(0, NULL);\
-        size_t m = (ol)->m, n = (ol)->n;\
-        for(size_t i=0, r=0; i < m; i++, r+=inc) {\
+        for(size_t i=0; i < m; i++) {\
             mp_obj_t row = mp_obj_new_list(n, NULL);\
             mp_obj_list_t *row_ptr = MP_OBJ_TO_PTR(row);\
-            for(size_t j=0, s=0; j < n; j++, s+=inc) {\
-                row_ptr->items[j] = mp_const_false;\
-                if((op) == MP_BINARY_OP_LESS) {\
-                    if(left[i*n+j] < right[r*n+s]) row_ptr->items[j] = mp_const_true;\
-                } else if((op) == MP_BINARY_OP_LESS_EQUAL) {\
-                    if(left[i*n+j] <= right[r*n+s]) row_ptr->items[j] = mp_const_true;\
-                } else if((op) == MP_BINARY_OP_MORE) {\
-                    if(left[i*n+j] > right[r*n+s]) row_ptr->items[j] = mp_const_true;\
-                } else if((op) == MP_BINARY_OP_MORE_EQUAL) {\
-                    if(left[i*n+j] >= right[r*n+s]) row_ptr->items[j] = mp_const_true;\
-                }\
-            }\
+			if((op) == MP_BINARY_OP_LESS) { for(size_t j=0; j < n; j++, left+=linc, right+=rinc) row_ptr->items[j] = *left < *right ? mp_const_true : mp_const_false; }\
+			else if((op) == MP_BINARY_OP_LESS_EQUAL) { for(size_t j=0; j < n; j++, left+=linc, right+=rinc) row_ptr->items[j] = *left <= *right ? mp_const_true : mp_const_false; }\
+			else if((op) == MP_BINARY_OP_MORE) { for(size_t j=0; j < n; j++, left+=linc, right+=rinc) row_ptr->items[j] = *left > *right ? mp_const_true : mp_const_false; }\
+			else if((op) == MP_BINARY_OP_MORE_EQUAL) { for(size_t j=0; j < n; j++, left+=linc, right+=rinc) row_ptr->items[j] = *left >= *right ? mp_const_true : mp_const_false; }\
+			else if((op) == MP_BINARY_OP_EQUAL) { for(size_t j=0; j < n; j++, left+=linc, right+=rinc) row_ptr->items[j] = *left == *right ? mp_const_true : mp_const_false; }\
+			else if((op) == MP_BINARY_OP_NOT_EQUAL) { for(size_t j=0; j < n; j++, left+=linc, right+=rinc) row_ptr->items[j] = *left != *right ? mp_const_true : mp_const_false; }\
             if(m == 1) return row;\
             mp_obj_list_append(out_list, row);\
-        }\
+		}\
         return out_list;\
     }\
 } while(0)

code/ndarray_properties.h

@@ -20,43 +20,42 @@
 
 #include "ndarray.h"
 
+#if CIRCUITPY
 typedef struct _mp_obj_property_t {
     mp_obj_base_t base;
     mp_obj_t proxy[3]; // getter, setter, deleter
 } mp_obj_property_t;
 
-/* v923z: it is not at all clear to me, why this must be declared; it should already be in obj.h */
-typedef struct _mp_obj_none_t {
-    mp_obj_base_t base;
-} mp_obj_none_t;
-
-const mp_obj_type_t mp_type_NoneType;
-const mp_obj_none_t mp_const_none_obj = {{&mp_type_NoneType}};
-
 MP_DEFINE_CONST_FUN_OBJ_1(ndarray_get_shape_obj, ndarray_shape);
 MP_DEFINE_CONST_FUN_OBJ_1(ndarray_get_size_obj, ndarray_size);
 MP_DEFINE_CONST_FUN_OBJ_1(ndarray_get_itemsize_obj, ndarray_itemsize);
-MP_DEFINE_CONST_FUN_OBJ_KW(ndarray_flatten_obj, 1, ndarray_flatten);
 
 STATIC const mp_obj_property_t ndarray_shape_obj = {
     .base.type = &mp_type_property,
     .proxy = {(mp_obj_t)&ndarray_get_shape_obj,
-              (mp_obj_t)&mp_const_none_obj,
-              (mp_obj_t)&mp_const_none_obj},
+              mp_const_none,
+              mp_const_none },
 };
 
 STATIC const mp_obj_property_t ndarray_size_obj = {
     .base.type = &mp_type_property,
     .proxy = {(mp_obj_t)&ndarray_get_size_obj,
-              (mp_obj_t)&mp_const_none_obj,
-              (mp_obj_t)&mp_const_none_obj},
+              mp_const_none,
+              mp_const_none },
 };
 
 STATIC const mp_obj_property_t ndarray_itemsize_obj = {
     .base.type = &mp_type_property,
     .proxy = {(mp_obj_t)&ndarray_get_itemsize_obj,
-              (mp_obj_t)&mp_const_none_obj,
-              (mp_obj_t)&mp_const_none_obj},
+              mp_const_none,
+              mp_const_none },
 };
+#else
+
+MP_DEFINE_CONST_FUN_OBJ_1(ndarray_size_obj, ndarray_size);
+MP_DEFINE_CONST_FUN_OBJ_1(ndarray_itemsize_obj, ndarray_itemsize);
+MP_DEFINE_CONST_FUN_OBJ_1(ndarray_shape_obj, ndarray_shape);
+
+#endif
 
 #endif

code/numerical.c

@@ -31,58 +31,7 @@ enum NUMERICAL_FUNCTION_TYPE {
     NUMERICAL_STD,
 };
 
-mp_obj_t numerical_linspace(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_const_none } },
-        { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
-        { MP_QSTR_num, MP_ARG_INT, {.u_int = 50} },
-        { MP_QSTR_endpoint, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_true} },
-        { MP_QSTR_retstep, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_false} },
-        { MP_QSTR_dtype, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = NDARRAY_FLOAT} },
-    };
-
-    mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
-
-    uint16_t len = args[2].u_int;
-    if(len < 2) {
-        mp_raise_ValueError(translate("number of points must be at least 2"));
-    }
-    mp_float_t value, step;
-    value = mp_obj_get_float(args[0].u_obj);
-    uint8_t typecode = args[5].u_int;
-    if(args[3].u_obj == mp_const_true) step = (mp_obj_get_float(args[1].u_obj)-value)/(len-1);
-    else step = (mp_obj_get_float(args[1].u_obj)-value)/len;
-    ndarray_obj_t *ndarray = create_new_ndarray(1, len, typecode);
-    if(typecode == NDARRAY_UINT8) {
-        uint8_t *array = (uint8_t *)ndarray->array->items;
-        for(size_t i=0; i < len; i++, value += step) array[i] = (uint8_t)value;
-    } else if(typecode == NDARRAY_INT8) {
-        int8_t *array = (int8_t *)ndarray->array->items;
-        for(size_t i=0; i < len; i++, value += step) array[i] = (int8_t)value;
-    } else if(typecode == NDARRAY_UINT16) {
-        uint16_t *array = (uint16_t *)ndarray->array->items;
-        for(size_t i=0; i < len; i++, value += step) array[i] = (uint16_t)value;
-    } else if(typecode == NDARRAY_INT16) {
-        int16_t *array = (int16_t *)ndarray->array->items;
-        for(size_t i=0; i < len; i++, value += step) array[i] = (int16_t)value;
-    } else {
-        mp_float_t *array = (mp_float_t *)ndarray->array->items;
-        for(size_t i=0; i < len; i++, value += step) array[i] = value;
-    }
-    if(args[4].u_obj == mp_const_false) {
-        return MP_OBJ_FROM_PTR(ndarray);
-    } else {
-        mp_obj_t tuple[2];
-        tuple[0] = ndarray;
-        tuple[1] = mp_obj_new_float(step);
-        return mp_obj_new_tuple(2, tuple);
-    }
-}
-
-MP_DEFINE_CONST_FUN_OBJ_KW(numerical_linspace_obj, 2, numerical_linspace);
-
-void axis_sorter(ndarray_obj_t *ndarray, mp_obj_t axis, size_t *m, size_t *n, size_t *N, 
+static void axis_sorter(ndarray_obj_t *ndarray, mp_obj_t axis, size_t *m, size_t *n, size_t *N, 
                  size_t *increment, size_t *len, size_t *start_inc) {
     if(axis == mp_const_none) { // flatten the array
         *m = 1;
@@ -108,7 +57,7 @@ void axis_sorter(ndarray_obj_t *ndarray, mp_obj_t axis, size_t *m, size_t *n, si
     }    
 }
 
-mp_obj_t numerical_sum_mean_std_iterable(mp_obj_t oin, uint8_t optype, size_t ddof) {
+static mp_obj_t numerical_sum_mean_std_iterable(mp_obj_t oin, uint8_t optype, size_t ddof) {
     mp_float_t value, sum = 0.0, sq_sum = 0.0;
     mp_obj_iter_buf_t iter_buf;
     mp_obj_t item, iterable = mp_getiter(oin, &iter_buf);
@@ -166,7 +115,7 @@ STATIC mp_obj_t numerical_sum_mean_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis
     return MP_OBJ_FROM_PTR(results);
 }
 
-mp_obj_t numerical_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, size_t ddof) {
+static mp_obj_t numerical_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, size_t ddof) {
     size_t m, n, increment, start, start_inc, N, len; 
     mp_float_t sum, sum_sq;
     
@@ -199,20 +148,32 @@ mp_obj_t numerical_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, size_t ddo
     return MP_OBJ_FROM_PTR(results);
 }
 
-mp_obj_t numerical_argmin_argmax_iterable(mp_obj_t oin, mp_obj_t axis, uint8_t optype) {
+static mp_obj_t numerical_argmin_argmax_iterable(mp_obj_t oin, uint8_t optype) {
+	if(MP_OBJ_SMALL_INT_VALUE(mp_obj_len_maybe(oin)) == 0) {
+		mp_raise_ValueError(translate("attempt to get argmin/argmax of an empty sequence"));
+	}
     size_t idx = 0, best_idx = 0;
     mp_obj_iter_buf_t iter_buf;
     mp_obj_t iterable = mp_getiter(oin, &iter_buf);
-    mp_obj_t best_obj = MP_OBJ_NULL;
     mp_obj_t item;
-    mp_uint_t op = MP_BINARY_OP_LESS;
-    if((optype == NUMERICAL_ARGMAX) || (optype == NUMERICAL_MAX)) op = MP_BINARY_OP_MORE;
-    while ((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
-        if ((best_obj == MP_OBJ_NULL) || (mp_binary_op(op, item, best_obj) == mp_const_true)) {
-            best_obj = item;
-            best_idx = idx;
-        }
-        idx++;
+    uint8_t op = 0; // argmin, min
+    if((optype == NUMERICAL_ARGMAX) || (optype == NUMERICAL_MAX)) op = 1;
+    item = mp_iternext(iterable);
+	mp_obj_t best_obj = item;
+    mp_float_t value, best_value = mp_obj_get_float(item);
+    value = best_value;
+    while((item = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
+		idx++;
+		value = mp_obj_get_float(item);
+		if((op == 0) && (value < best_value)) {
+			best_obj = item;
+			best_idx = idx;
+			best_value = value;
+		} else if((op == 1) && (value > best_value)) {
+			best_obj = item;
+			best_idx = idx;
+			best_value = value;
+		}
     }
     if((optype == NUMERICAL_ARGMIN) || (optype == NUMERICAL_ARGMAX)) {
         return MP_OBJ_NEW_SMALL_INT(best_idx);
@@ -221,7 +182,7 @@ mp_obj_t numerical_argmin_argmax_iterable(mp_obj_t oin, mp_obj_t axis, uint8_t o
     }    
 }
 
-mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {
+static mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype) {
     size_t m, n, increment, start, start_inc, N, len;
     axis_sorter(ndarray, axis, &m, &n, &N, &increment, &len, &start_inc);
     ndarray_obj_t *results;
@@ -229,6 +190,9 @@ mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis,
         // we could save some RAM by taking NDARRAY_UINT8, if the dimensions 
         // are smaller than 256, but the code would become more involving 
         // (we would also need extra flash space)
+        if(ndarray->array->len == 0) {
+			mp_raise_ValueError(translate("attempt to get argmin/argmax of an empty sequence"));
+		}
         results = create_new_ndarray(m, n, NDARRAY_UINT16);
     } else {
         results = create_new_ndarray(m, n, ndarray->array->typecode);
@@ -236,22 +200,35 @@ mp_obj_t numerical_argmin_argmax_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis,
     
     for(size_t j=0; j < N; j++) { // result index
         start = j * start_inc;
-        if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {
-            if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {
+        if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {
+            if((ndarray->array->typecode == NDARRAY_UINT8)) {
                 RUN_ARGMIN(ndarray, results, uint8_t, uint8_t, len, start, increment, optype, j);
+            } else if((ndarray->array->typecode == NDARRAY_INT8)) {
+                RUN_ARGMIN(ndarray, results, int8_t, int8_t, len, start, increment, optype, j);
+            } if((ndarray->array->typecode == NDARRAY_UINT16)) {
+                RUN_ARGMIN(ndarray, results, uint16_t, uint16_t, len, start, increment, optype, j);
+            } else if((ndarray->array->typecode == NDARRAY_INT16)) {
+                RUN_ARGMIN(ndarray, results, int16_t, int16_t, len, start, increment, optype, j);
             } else {
-                RUN_ARGMIN(ndarray, results, uint8_t, uint16_t, len, start, increment, optype, j);                
-            }
-        } else if((ndarray->array->typecode == NDARRAY_UINT16) || (ndarray->array->typecode == NDARRAY_INT16)) {
-            RUN_ARGMIN(ndarray, results, uint16_t, uint16_t, len, start, increment, optype, j);
-        } else {
-            if((optype == NUMERICAL_MAX) || (optype == NUMERICAL_MIN)) {
                 RUN_ARGMIN(ndarray, results, mp_float_t, mp_float_t, len, start, increment, optype, j);
+            }
+        } else {
+            if((ndarray->array->typecode == NDARRAY_UINT8)) {
+                RUN_ARGMIN(ndarray, results, uint8_t, uint16_t, len, start, increment, optype, j);
+            } else if((ndarray->array->typecode == NDARRAY_INT8)) {
+                RUN_ARGMIN(ndarray, results, int8_t, uint16_t, len, start, increment, optype, j);
+            } if((ndarray->array->typecode == NDARRAY_UINT16)) {
+                RUN_ARGMIN(ndarray, results, uint16_t, uint16_t, len, start, increment, optype, j);
+            } else if((ndarray->array->typecode == NDARRAY_INT16)) {
+                RUN_ARGMIN(ndarray, results, int16_t, uint16_t, len, start, increment, optype, j);
             } else {
-                RUN_ARGMIN(ndarray, results, mp_float_t, uint16_t, len, start, increment, optype, j);                
+                RUN_ARGMIN(ndarray, results, mp_float_t, uint16_t, len, start, increment, optype, j);
             }
         }
     }
+	if(results->array->len == 1) {
+		return mp_binary_get_val_array(results->array->typecode, results->array->items, 0);
+	}
     return MP_OBJ_FROM_PTR(results);
 }
 
@@ -278,7 +255,7 @@ STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_m
             case NUMERICAL_ARGMIN:
             case NUMERICAL_MAX:
             case NUMERICAL_ARGMAX:
-                return numerical_argmin_argmax_iterable(oin, axis, optype);
+                return numerical_argmin_argmax_iterable(oin, optype);
             case NUMERICAL_SUM:
             case NUMERICAL_MEAN:
                 return numerical_sum_mean_std_iterable(oin, optype, 0);
@@ -305,43 +282,43 @@ STATIC mp_obj_t numerical_function(size_t n_args, const mp_obj_t *pos_args, mp_m
     return mp_const_none;
 }
 
-mp_obj_t numerical_min(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_min(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MIN);
 }
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_min_obj, 1, numerical_min);
 
-mp_obj_t numerical_max(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_max(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MAX);
 }
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_max_obj, 1, numerical_max);
 
-mp_obj_t numerical_argmin(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_argmin(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     return numerical_function(n_args, pos_args, kw_args, NUMERICAL_ARGMIN);
 }
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argmin_obj, 1, numerical_argmin);
 
-mp_obj_t numerical_argmax(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_argmax(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     return numerical_function(n_args, pos_args, kw_args, NUMERICAL_ARGMAX);
 }
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argmax_obj, 1, numerical_argmax);
 
-mp_obj_t numerical_sum(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_sum(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     return numerical_function(n_args, pos_args, kw_args, NUMERICAL_SUM);
 }
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sum_obj, 1, numerical_sum);
 
-mp_obj_t numerical_mean(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_mean(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
     return numerical_function(n_args, pos_args, kw_args, NUMERICAL_MEAN);
 }
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_mean_obj, 1, numerical_mean);
 
-mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_std(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_const_none } } ,
         { MP_QSTR_axis, MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
@@ -371,7 +348,7 @@ mp_obj_t numerical_std(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_arg
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_std_obj, 1, numerical_std);
 
-mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_roll(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_const_none  } },
         { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
@@ -379,7 +356,7 @@ mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(2, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
     
     mp_obj_t oin = args[0].u_obj;
     int16_t shift = mp_obj_get_int(args[1].u_obj);
@@ -453,14 +430,14 @@ mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_roll_obj, 2, numerical_roll);
 
-mp_obj_t numerical_flip(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_flip(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_const_none } },
         { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none } },
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
     
     if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {
         mp_raise_TypeError(translate("flip argument must be an ndarray"));
@@ -502,7 +479,7 @@ mp_obj_t numerical_flip(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_flip_obj, 1, numerical_flip);
 
-mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_diff(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_const_none } },
         { MP_QSTR_n, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 1 } },
@@ -510,7 +487,7 @@ mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
     
     if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {
         mp_raise_TypeError(translate("diff argument must be an ndarray"));
@@ -572,7 +549,7 @@ mp_obj_t numerical_diff(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_diff_obj, 1, numerical_diff);
 
-mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inplace) {
+static mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inplace) {
     if(!MP_OBJ_IS_TYPE(oin, &ulab_ndarray_type)) {
         mp_raise_TypeError(translate("sort argument must be an ndarray"));
     }
@@ -602,7 +579,7 @@ mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inplace) {
     } else if(mp_obj_get_int(axis) == 0) { // sort along vertical axis
         increment = ndarray->n;
         start_inc = 1;
-        end = ndarray->m;
+        end = ndarray->n;
         N = ndarray->m;
     } else {
         mp_raise_ValueError(translate("axis must be -1, 0, None, or 1"));
@@ -611,7 +588,7 @@ mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inplace) {
     size_t q, k, p, c;
 
     for(size_t start=0; start < end; start+=start_inc) {
-        q = N; 
+        q = N;
         k = (q >> 1);
         if((ndarray->array->typecode == NDARRAY_UINT8) || (ndarray->array->typecode == NDARRAY_INT8)) {
             HEAPSORT(uint8_t, ndarray);
@@ -629,14 +606,14 @@ mp_obj_t numerical_sort_helper(mp_obj_t oin, mp_obj_t axis, uint8_t inplace) {
 }
 
 // numpy function
-mp_obj_t numerical_sort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_sort(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_const_none } },
         { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_int = -1 } },
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
 
     return numerical_sort_helper(args[0].u_obj, args[1].u_obj, 0);
 }
@@ -644,27 +621,27 @@ mp_obj_t numerical_sort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sort_obj, 1, numerical_sort);
 
 // method of an ndarray
-mp_obj_t numerical_sort_inplace(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_sort_inplace(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_const_none } },
         { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_int = -1 } },
     };
 
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
 
     return numerical_sort_helper(args[0].u_obj, args[1].u_obj, 1);
 }
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_sort_inplace_obj, 1, numerical_sort_inplace);
 
-mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) {
+static mp_obj_t numerical_argsort(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_const_none } },
         { MP_QSTR_axis, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_int = -1 } },
     };
     mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)];
-    mp_arg_parse_all(1, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
+    mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args);
     if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {
         mp_raise_TypeError(translate("argsort argument must be an ndarray"));
     }
@@ -674,8 +651,6 @@ mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
     if(args[1].u_obj == mp_const_none) { // flatten the array
         m = 1;
         n = ndarray->array->len;
-        ndarray->m = m;
-        ndarray->n = n;
         increment = 1;
         start_inc = ndarray->n;
         end = ndarray->n;
@@ -693,7 +668,7 @@ mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
         n = ndarray->n;
         increment = n;
         start_inc = 1;
-        end = m;
+        end = n;
         N = m;
     } else {
         mp_raise_ValueError(translate("axis must be -1, 0, None, or 1"));
@@ -730,9 +705,8 @@ mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
 
 MP_DEFINE_CONST_FUN_OBJ_KW(numerical_argsort_obj, 1, numerical_argsort);
 
-#if !CIRCUITPY
 STATIC const mp_rom_map_elem_t ulab_numerical_globals_table[] = {
-    { MP_OBJ_NEW_QSTR(MP_QSTR_linspace), (mp_obj_t)&numerical_linspace_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_numerical) },
     { MP_OBJ_NEW_QSTR(MP_QSTR_sum), (mp_obj_t)&numerical_sum_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_mean), (mp_obj_t)&numerical_mean_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_std), (mp_obj_t)&numerical_std_obj },
@@ -753,6 +727,5 @@ mp_obj_module_t ulab_numerical_module = {
     .base = { &mp_type_module },
     .globals = (mp_obj_dict_t*)&mp_module_ulab_numerical_globals,
 };
-#endif
 
 #endif

code/numerical.h

@@ -19,9 +19,7 @@
 
 extern mp_obj_module_t ulab_numerical_module;
 
-// TODO: implement minimum/maximum, and cumsum
-//mp_obj_t numerical_minimum(mp_obj_t , mp_obj_t );
-//mp_obj_t numerical_maximum(mp_obj_t , mp_obj_t );
+// TODO: implement cumsum
 //mp_obj_t numerical_cumsum(size_t , const mp_obj_t *, mp_map_t *);
 
 #define RUN_ARGMIN(in, out, typein, typeout, len, start, increment, op, pos) do {\
@@ -32,13 +30,13 @@ extern mp_obj_module_t ulab_numerical_module;
         for(size_t i=1; i < (len); i++) {\
             if(array[(start)+i*(increment)] > array[(start)+best_index*(increment)]) best_index = i;\
         }\
-        if((op) == NUMERICAL_MAX) outarray[(pos)] = array[(start)+best_index*(increment)];\
+        if((op) == NUMERICAL_MAX) outarray[(pos)] = (typeout)array[(start)+best_index*(increment)];\
         else outarray[(pos)] = best_index;\
     } else{\
         for(size_t i=1; i < (len); i++) {\
             if(array[(start)+i*(increment)] < array[(start)+best_index*(increment)]) best_index = i;\
         }\
-        if((op) == NUMERICAL_MIN) outarray[(pos)] = array[(start)+best_index*(increment)];\
+        if((op) == NUMERICAL_MIN) outarray[(pos)] = (typeout)array[(start)+best_index*(increment)];\
         else outarray[(pos)] = best_index;\
     }\
 } while(0)
@@ -46,8 +44,8 @@ extern mp_obj_module_t ulab_numerical_module;
 #define RUN_SUM(ndarray, type, optype, len, start, increment) do {\
     type *array = (type *)(ndarray)->array->items;\
     type value;\
-    for(size_t j=0; j < (len); j++) {\
-        value = array[(start)+j*(increment)];\
+    for(size_t k=0; k < (len); k++) {\
+        value = array[(start)+k*(increment)];\
         sum += value;\
     }\
 } while(0)
@@ -55,12 +53,12 @@ extern mp_obj_module_t ulab_numerical_module;
 #define RUN_STD(ndarray, type, len, start, increment) do {\
     type *array = (type *)(ndarray)->array->items;\
     mp_float_t value;\
-    for(size_t j=0; j < (len); j++) {\
-        sum += array[(start)+j*(increment)];\
+    for(size_t k=0; k < (len); k++) {\
+        sum += array[(start)+k*(increment)];\
     }\
     sum /= (len);\
-    for(size_t j=0; j < (len); j++) {\
-        value = (array[(start)+j*(increment)] - sum);\
+    for(size_t k=0; k < (len); k++) {\
+        value = (array[(start)+k*(increment)] - sum);\
         sum_sq += value * value;\
     }\
 } while(0)
@@ -79,33 +77,34 @@ extern mp_obj_module_t ulab_numerical_module;
 
 #define HEAPSORT(type, ndarray) do {\
     type *array = (type *)(ndarray)->array->items;\
+    array += start;\
     type tmp;\
     for (;;) {\
         if (k > 0) {\
-            tmp = array[start+(--k)*increment];\
+            tmp = array[(--k)*increment];\
         } else {\
             q--;\
             if(q == 0) {\
                 break;\
             }\
-            tmp = array[start+q*increment];\
-            array[start+q*increment] = array[start];\
+            tmp = array[q*increment];\
+            array[q*increment] = array[0];\
         }\
         p = k;\
         c = k + k + 1;\
         while (c < q) {\
-            if((c + 1 < q)  &&  (array[start+(c+1)*increment] > array[start+c*increment])) {\
+            if((c + 1 < q)  &&  (array[(c+1)*increment] > array[c*increment])) {\
                 c++;\
             }\
-            if(array[start+c*increment] > tmp) {\
-                array[start+p*increment] = array[start+c*increment];\
+            if(array[c*increment] > tmp) {\
+                array[p*increment] = array[c*increment];\
                 p = c;\
                 c = p + p + 1;\
             } else {\
                 break;\
             }\
         }\
-        array[start+p*increment] = tmp;\
+        array[p*increment] = tmp;\
     }\
 } while(0)
 
@@ -148,7 +147,6 @@ extern mp_obj_module_t ulab_numerical_module;
     }\
 } while(0)
 
-MP_DECLARE_CONST_FUN_OBJ_KW(numerical_linspace_obj);
 MP_DECLARE_CONST_FUN_OBJ_KW(numerical_min_obj);
 MP_DECLARE_CONST_FUN_OBJ_KW(numerical_max_obj);
 MP_DECLARE_CONST_FUN_OBJ_KW(numerical_argmin_obj);

code/poly.c

@@ -6,7 +6,9 @@
  *
  * The MIT License (MIT)
  *
- * Copyright (c) 2019-2020 Zoltán Vörös
+ * Copyright (c) 2020 Jeff Epler for Adafruit Industries
+ *               2019-2020 Zoltán Vörös
+ *
 */
 
 #include "py/obj.h"
@@ -17,26 +19,7 @@
 #include "poly.h"
 
 #if ULAB_POLY_MODULE
-bool object_is_nditerable(mp_obj_t o_in) {
-    if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type) || 
-      MP_OBJ_IS_TYPE(o_in, &mp_type_tuple) || 
-      MP_OBJ_IS_TYPE(o_in, &mp_type_list) || 
-      MP_OBJ_IS_TYPE(o_in, &mp_type_range)) {
-        return true;
-    }
-    return false;
-}
-
-size_t get_nditerable_len(mp_obj_t o_in) {
-    if(MP_OBJ_IS_TYPE(o_in, &ulab_ndarray_type)) {
-        ndarray_obj_t *in = MP_OBJ_TO_PTR(o_in);
-        return in->array->len;
-    } else {
-        return (size_t)mp_obj_get_int(mp_obj_len_maybe(o_in));
-    }
-}
-
-mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
+static mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
     // TODO: return immediately, if o_p is not an iterable
     // TODO: there is a bug here: matrices won't work, 
     // because there is a single iteration loop
@@ -85,11 +68,11 @@ mp_obj_t poly_polyval(mp_obj_t o_p, mp_obj_t o_x) {
 
 MP_DEFINE_CONST_FUN_OBJ_2(poly_polyval_obj, poly_polyval);
 
-mp_obj_t poly_polyfit(size_t  n_args, const mp_obj_t *args) {
+static mp_obj_t poly_polyfit(size_t  n_args, const mp_obj_t *args) {
     if((n_args != 2) && (n_args != 3)) {
         mp_raise_ValueError(translate("number of arguments must be 2, or 3"));
     }
-    if(!object_is_nditerable(args[0])) {
+    if(!ndarray_object_is_nditerable(args[0])) {
         mp_raise_ValueError(translate("input data must be an iterable"));
     }
     uint16_t lenx = 0, leny = 0;
@@ -111,9 +94,12 @@ mp_obj_t poly_polyfit(size_t  n_args, const mp_obj_t *args) {
         }
         y = m_new(mp_float_t, leny);
         fill_array_iterable(y, args[0]);
-    } else if(n_args == 3) {
+    } else /* n_args == 3 */ {
+        if(!ndarray_object_is_nditerable(args[1])) {
+            mp_raise_ValueError(translate("input data must be an iterable"));
+        }
         lenx = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));
-        leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));
+        leny = (uint16_t)mp_obj_get_int(mp_obj_len_maybe(args[1]));
         if(lenx != leny) {
             mp_raise_ValueError(translate("input vectors must be of equal length"));
         }
@@ -197,7 +183,6 @@ mp_obj_t poly_polyfit(size_t  n_args, const mp_obj_t *args) {
 
 MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(poly_polyfit_obj, 2, 3, poly_polyfit);
 
-#if !CIRCUITPY
 STATIC const mp_rom_map_elem_t ulab_poly_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_poly) },
     { MP_OBJ_NEW_QSTR(MP_QSTR_polyval), (mp_obj_t)&poly_polyval_obj },
@@ -210,6 +195,5 @@ mp_obj_module_t ulab_poly_module = {
     .base = { &mp_type_module },
     .globals = (mp_obj_dict_t*)&mp_module_ulab_poly_globals,
 };
-#endif
 
 #endif

code/ulab.c

@@ -21,30 +21,40 @@
 #include "ulab.h"
 #include "ndarray.h"
 #include "ndarray_properties.h"
+#include "create.h"
 #include "linalg.h"
 #include "vectorise.h"
 #include "poly.h"
 #include "fft.h"
 #include "filter.h"
 #include "numerical.h"
+#include "compare.h"
+#include "approx.h"
 #include "extras.h"
 
-STATIC MP_DEFINE_STR_OBJ(ulab_version_obj, "0.34.0");
+STATIC MP_DEFINE_STR_OBJ(ulab_version_obj, "0.50.1");
+
+MP_DEFINE_CONST_FUN_OBJ_KW(ndarray_flatten_obj, 1, ndarray_flatten);
 
 STATIC const mp_rom_map_elem_t ulab_ndarray_locals_dict_table[] = {
-    { MP_ROM_QSTR(MP_QSTR_flatten), MP_ROM_PTR(&ndarray_flatten_obj) },
     { MP_ROM_QSTR(MP_QSTR_reshape), MP_ROM_PTR(&ndarray_reshape_obj) },
     { MP_ROM_QSTR(MP_QSTR_transpose), MP_ROM_PTR(&ndarray_transpose_obj) },
+    { MP_ROM_QSTR(MP_QSTR_flatten), MP_ROM_PTR(&ndarray_flatten_obj) },
     { MP_ROM_QSTR(MP_QSTR_shape), MP_ROM_PTR(&ndarray_shape_obj) },
     { MP_ROM_QSTR(MP_QSTR_size), MP_ROM_PTR(&ndarray_size_obj) },
     { MP_ROM_QSTR(MP_QSTR_itemsize), MP_ROM_PTR(&ndarray_itemsize_obj) },
-//    { MP_ROM_QSTR(MP_QSTR_sort), MP_ROM_PTR(&numerical_sort_inplace_obj) },
+#if CIRCUITPY
+    { MP_ROM_QSTR(MP_QSTR_sort), MP_ROM_PTR(&numerical_sort_inplace_obj) },
+#endif
 };
 
 STATIC MP_DEFINE_CONST_DICT(ulab_ndarray_locals_dict, ulab_ndarray_locals_dict_table);
 
 const mp_obj_type_t ulab_ndarray_type = {
     { &mp_type_type },
+#if defined(MP_TYPE_FLAG_EQ_CHECKS_OTHER_TYPE) && defined(MP_TYPE_FLAG_EQ_HAS_NEQ_TEST)
+    .flags = MP_TYPE_FLAG_EQ_CHECKS_OTHER_TYPE | MP_TYPE_FLAG_EQ_HAS_NEQ_TEST,
+#endif
     .name = MP_QSTR_ndarray,
     .print = ndarray_print,
     .make_new = ndarray_make_new,
@@ -56,11 +66,16 @@ const mp_obj_type_t ulab_ndarray_type = {
     .locals_dict = (mp_obj_dict_t*)&ulab_ndarray_locals_dict,
 };
 
-#if !CIRCUITPY
 STATIC const mp_map_elem_t ulab_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_ulab) },
     { MP_ROM_QSTR(MP_QSTR___version__), MP_ROM_PTR(&ulab_version_obj) },
+    { MP_ROM_QSTR(MP_QSTR_set_printoptions), (mp_obj_t)&ndarray_set_printoptions_obj },
+    { MP_ROM_QSTR(MP_QSTR_get_printoptions), (mp_obj_t)&ndarray_get_printoptions_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_array), (mp_obj_t)&ulab_ndarray_type },
+    { MP_ROM_QSTR(MP_QSTR_zeros), (mp_obj_t)&create_zeros_obj },
+    { MP_ROM_QSTR(MP_QSTR_ones), (mp_obj_t)&create_ones_obj },
+    { MP_ROM_QSTR(MP_QSTR_eye), (mp_obj_t)&create_eye_obj },
+    { MP_ROM_QSTR(MP_QSTR_linspace), (mp_obj_t)&create_linspace_obj },
     #if ULAB_LINALG_MODULE
     { MP_ROM_QSTR(MP_QSTR_linalg), MP_ROM_PTR(&ulab_linalg_module) },
     #endif
@@ -79,6 +94,12 @@ STATIC const mp_map_elem_t ulab_globals_table[] = {
     #if ULAB_FILTER_MODULE
     { MP_ROM_QSTR(MP_QSTR_filter), MP_ROM_PTR(&ulab_filter_module) },
     #endif
+    #if ULAB_COMPARE_MODULE
+    { MP_ROM_QSTR(MP_QSTR_compare), MP_ROM_PTR(&ulab_compare_module) },
+    #endif
+    #if ULAB_APPROX_MODULE
+    { MP_ROM_QSTR(MP_QSTR_approx), MP_ROM_PTR(&ulab_approx_module) },
+    #endif
     #if ULAB_EXTRAS_MODULE
     { MP_ROM_QSTR(MP_QSTR_extras), MP_ROM_PTR(&ulab_extras_module) },
     #endif
@@ -101,4 +122,3 @@ mp_obj_module_t ulab_user_cmodule = {
 };
 
 MP_REGISTER_MODULE(MP_QSTR_ulab, ulab_user_cmodule, MODULE_ULAB_ENABLED);
-#endif

code/ulab.h

@@ -12,25 +12,52 @@
 #ifndef __ULAB__
 #define __ULAB__
 
-// vectorise (all functions) takes approx. 3 kB of flash space
+// the create module is always included
+#define ULAB_CREATE_MODULE (1)
+
+// vectorise (all functions) takes approx. 6 kB of flash space
+#ifndef ULAB_VECTORISE_MODULE
 #define ULAB_VECTORISE_MODULE (1)
+#endif
 
-// linalg adds around 6 kB
+// linalg adds around 6 kB to the firmware
+#ifndef ULAB_LINALG_MODULE
 #define ULAB_LINALG_MODULE (1)
+#endif
 
-// poly is approx. 2.5 kB
+// poly requires approx. 2.5 kB
+#ifndef ULAB_POLY_MODULE
 #define ULAB_POLY_MODULE (1)
+#endif
 
 // numerical is about 12 kB
+#ifndef ULAB_NUMERICAL_MODULE
 #define ULAB_NUMERICAL_MODULE (1)
+#endif
 
 // FFT costs about 2 kB of flash space
+#ifndef ULAB_FFT_MODULE
 #define ULAB_FFT_MODULE (1)
+#endif
 
-// the filter module takes about 1 kB of flash space
+// the filter module occupies about 1 kB of flash space
+#ifndef ULAB_FILTER_MODULE
 #define ULAB_FILTER_MODULE (1)
+#endif
+
+// the compare module consumes about 4 kB of flash space
+#ifndef ULAB_COMPARE_MODULE
+#define ULAB_COMPARE_MODULE (1)
+#endif
+
+// the approx module consumes about ... kB of flash space
+#ifndef ULAB_APPROX_MODULE
+#define ULAB_APPROX_MODULE (1)
+#endif
 
 // user-defined modules
-#define ULAB_EXTRAS_MODULE (0)
+#ifndef ULAB_EXTRAS_MODULE
+#define ULAB_EXTRAS_MODULE (1)
+#endif
 
 #endif

code/vectorise.c

@@ -23,7 +23,7 @@
 #endif
     
 #if ULAB_VECTORISE_MODULE
-mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {
+static mp_obj_t vectorise_generic_vector(mp_obj_t o_in, mp_float_t (*f)(mp_float_t)) {
     // Return a single value, if o_in is not iterable
     if(mp_obj_is_float(o_in) || MP_OBJ_IS_INT(o_in)) {
             return mp_obj_new_float(f(mp_obj_get_float(o_in)));
@@ -135,11 +135,171 @@ MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tan_obj, vectorise_tan);
 MATH_FUN_1(tanh, tanh);
 MP_DEFINE_CONST_FUN_OBJ_1(vectorise_tanh_obj, vectorise_tanh);
 
-#if !CIRCUITPY
+static mp_obj_t vectorise_around(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_const_none} },
+        { MP_QSTR_decimals, MP_ARG_KW_ONLY | MP_ARG_INT, {.u_int = 0 } }
+    };
+
+    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);
+	if(!MP_OBJ_IS_TYPE(args[0].u_obj, &ulab_ndarray_type)) {
+		mp_raise_TypeError(translate("first argument must be an ndarray"));
+	}
+    int8_t n = args[1].u_int;
+	mp_float_t mul = MICROPY_FLOAT_C_FUN(pow)(10.0, n);
+    ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
+    ndarray_obj_t *result = create_new_ndarray(ndarray->m, ndarray->n, NDARRAY_FLOAT);
+    mp_float_t *array = (mp_float_t *)result->array->items;
+    for(size_t i=0; i < ndarray->array->len; i++) {
+		mp_float_t f = ndarray_get_float_value(ndarray->array->items, ndarray->array->typecode, i);
+		*array++ = MICROPY_FLOAT_C_FUN(round)(f * mul) / mul;
+	}
+    return MP_OBJ_FROM_PTR(result);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(vectorise_around_obj, 1, vectorise_around);
+
+static mp_obj_t vectorise_arctan2(mp_obj_t x, mp_obj_t y) {
+	// the function is implemented for scalars and ndarrays only, with partial 
+	// broadcasting: arguments must be either scalars, or ndarrays of equal size/shape
+	if(!(MP_OBJ_IS_INT(x) || mp_obj_is_float(x) || MP_OBJ_IS_TYPE(x, &ulab_ndarray_type)) &&
+		!(MP_OBJ_IS_INT(y) || mp_obj_is_float(y) || MP_OBJ_IS_TYPE(y, &ulab_ndarray_type))) {
+		mp_raise_TypeError(translate("arctan2 is implemented for scalars and ndarrays only"));
+	}
+	ndarray_obj_t *ndarray_x, *ndarray_y;
+	if(MP_OBJ_IS_INT(x) || mp_obj_is_float(x)) {
+		ndarray_x = create_new_ndarray(1, 1, NDARRAY_FLOAT);
+		mp_float_t *array_x = (mp_float_t *)ndarray_x->array->items;
+		*array_x = mp_obj_get_float(x);
+	} else {
+		ndarray_x = MP_OBJ_TO_PTR(x);
+	}
+	if(MP_OBJ_IS_INT(y) || mp_obj_is_float(y)) {
+		ndarray_y = create_new_ndarray(1, 1, NDARRAY_FLOAT);
+		mp_float_t *array_y = (mp_float_t *)ndarray_y->array->items;
+		*array_y = mp_obj_get_float(y);
+	} else {
+		ndarray_y = MP_OBJ_TO_PTR(y);
+	}
+	// check, whether partial broadcasting is possible here
+	if((ndarray_x->m != ndarray_y->m) || (ndarray_x->n != ndarray_y->n)) {
+		if((ndarray_x->array->len != 1) && (ndarray_y->array->len != 1)) {
+            mp_raise_ValueError(translate("operands could not be broadcast together"));
+		}
+	}
+	size_t xinc = 0, yinc = 0;
+	size_t m = MAX(ndarray_x->m, ndarray_y->m);
+	size_t n = MAX(ndarray_x->n, ndarray_y->n);
+	size_t len = MAX(ndarray_x->array->len, ndarray_y->array->len);
+	if(ndarray_x->array->len != 1) {
+		xinc = 1;
+	}
+	if(ndarray_y->array->len != 1) {
+		yinc = 1;
+	}
+	size_t posx = 0, posy = 0;
+	ndarray_obj_t *result = create_new_ndarray(m, n, NDARRAY_FLOAT);
+	mp_float_t *array_r = (mp_float_t *)result->array->items;
+	for(size_t i=0; i < len; i++) {
+		mp_float_t value_x = ndarray_get_float_value(ndarray_x->array->items, ndarray_x->array->typecode, posx);
+		mp_float_t value_y = ndarray_get_float_value(ndarray_y->array->items, ndarray_y->array->typecode, posy);
+		*array_r++ = MICROPY_FLOAT_C_FUN(atan2)(value_x, value_y);
+		posx += xinc;
+		posy += yinc;
+	}
+    return MP_OBJ_FROM_PTR(result);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_2(vectorise_arctan2_obj, vectorise_arctan2);
+
+static mp_obj_t vectorise_vectorized_function_call(mp_obj_t self_in, size_t n_args, size_t n_kw, const mp_obj_t *args) {
+    (void) n_args;
+    (void) n_kw;
+    vectorized_function_obj_t *self = MP_OBJ_TO_PTR(self_in);
+    mp_obj_t avalue[1];
+    mp_obj_t fvalue;
+    if(MP_OBJ_IS_TYPE(args[0], &ulab_ndarray_type)) {
+        ndarray_obj_t *source = MP_OBJ_TO_PTR(args[0]);
+        ndarray_obj_t *target = create_new_ndarray(source->m, source->n, self->otypes);
+        for(size_t i=0; i < source->array->len; i++) {
+            avalue[0] = mp_binary_get_val_array(source->array->typecode, source->array->items, i);
+            fvalue = self->type->call(self->fun, 1, 0, avalue);
+            mp_binary_set_val_array(self->otypes, target->array->items, i, fvalue);
+        }
+        return MP_OBJ_FROM_PTR(target);
+    } else if(MP_OBJ_IS_TYPE(args[0], &mp_type_tuple) || MP_OBJ_IS_TYPE(args[0], &mp_type_list) ||
+        MP_OBJ_IS_TYPE(args[0], &mp_type_range)) { // i.e., the input is a generic iterable
+        size_t len = (size_t)mp_obj_get_int(mp_obj_len_maybe(args[0]));
+        ndarray_obj_t *target = create_new_ndarray(1, len, self->otypes);
+        mp_obj_iter_buf_t iter_buf;
+        mp_obj_t iterable = mp_getiter(args[0], &iter_buf);
+        size_t i=0;
+        while ((avalue[0] = mp_iternext(iterable)) != MP_OBJ_STOP_ITERATION) {
+            fvalue = self->type->call(self->fun, 1, 0, avalue);
+            mp_binary_set_val_array(self->otypes, target->array->items, i, fvalue);
+            i++;
+        }
+        return MP_OBJ_FROM_PTR(target);
+    } else if(mp_obj_is_int(args[0]) || mp_obj_is_float(args[0])) {
+        ndarray_obj_t *target = create_new_ndarray(1, 1, self->otypes);
+        fvalue = self->type->call(self->fun, 1, 0, args);
+        mp_binary_set_val_array(self->otypes, target->array->items, 0, fvalue);
+        return MP_OBJ_FROM_PTR(target);
+    } else {
+        mp_raise_ValueError(translate("wrong input type"));
+    }
+    return mp_const_none;
+}
+
+const mp_obj_type_t vectorise_function_type = {
+    { &mp_type_type },
+    .name = MP_QSTR_,
+    .call = vectorise_vectorized_function_call,
+};
+
+static mp_obj_t vectorise_vectorize(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_const_none} },
+        { MP_QSTR_otypes, MP_ARG_KW_ONLY | MP_ARG_OBJ, {.u_rom_obj = mp_const_none} }
+    };
+    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);
+    const mp_obj_type_t *type = mp_obj_get_type(args[0].u_obj);
+    if(type->call == NULL) {
+        mp_raise_TypeError(translate("first argument must be a callable"));
+    }
+    mp_obj_t _otypes = args[1].u_obj;
+    uint8_t otypes = NDARRAY_FLOAT;
+    if(_otypes == mp_const_none) {
+        // TODO: is this what numpy does?
+        otypes = NDARRAY_FLOAT;
+    } else if(mp_obj_is_int(_otypes)) {
+        otypes = mp_obj_get_int(_otypes);
+        if(otypes != NDARRAY_FLOAT && otypes != NDARRAY_UINT8 && otypes != NDARRAY_INT8 &&
+            otypes != NDARRAY_UINT16 && otypes != NDARRAY_INT16) {
+                mp_raise_ValueError(translate("wrong output type"));
+        }
+    }
+    else {
+        mp_raise_ValueError(translate("wrong output type"));
+    }
+    vectorized_function_obj_t *function = m_new_obj(vectorized_function_obj_t);
+    function->base.type = &vectorise_function_type;
+    function->otypes = otypes;
+    function->fun = args[0].u_obj;
+    function->type = type;
+    return MP_OBJ_FROM_PTR(function);
+}
+
+MP_DEFINE_CONST_FUN_OBJ_KW(vectorise_vectorize_obj, 1, vectorise_vectorize);
+
 STATIC const mp_rom_map_elem_t ulab_vectorise_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_vector) },
     { MP_OBJ_NEW_QSTR(MP_QSTR_acos), (mp_obj_t)&vectorise_acos_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_acosh), (mp_obj_t)&vectorise_acosh_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_arctan2), (mp_obj_t)&vectorise_arctan2_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_around), (mp_obj_t)&vectorise_around_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_asin), (mp_obj_t)&vectorise_asin_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_asinh), (mp_obj_t)&vectorise_asinh_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_atan), (mp_obj_t)&vectorise_atan_obj },
@@ -161,6 +321,7 @@ STATIC const mp_rom_map_elem_t ulab_vectorise_globals_table[] = {
     { MP_OBJ_NEW_QSTR(MP_QSTR_sqrt), (mp_obj_t)&vectorise_sqrt_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_tan), (mp_obj_t)&vectorise_tan_obj },
     { MP_OBJ_NEW_QSTR(MP_QSTR_tanh), (mp_obj_t)&vectorise_tanh_obj },
+    { MP_OBJ_NEW_QSTR(MP_QSTR_vectorize), (mp_obj_t)&vectorise_vectorize_obj },
 };
 
 STATIC MP_DEFINE_CONST_DICT(mp_module_ulab_vectorise_globals, ulab_vectorise_globals_table);
@@ -169,6 +330,5 @@ mp_obj_module_t ulab_vectorise_module = {
     .base = { &mp_type_module },
     .globals = (mp_obj_dict_t*)&mp_module_ulab_vectorise_globals,
 };
-#endif
 
 #endif

code/vectorise.h

@@ -17,6 +17,13 @@
 
 #if ULAB_VECTORISE_MODULE
 
+typedef struct _vectorized_function_obj_t {
+    mp_obj_base_t base;
+    uint8_t otypes;
+    mp_obj_t fun;
+    const mp_obj_type_t *type;
+} vectorized_function_obj_t;
+
 mp_obj_module_t ulab_vectorise_module;
 
 #define ITERATE_VECTOR(type, source, out) do {\
@@ -27,7 +34,7 @@ mp_obj_module_t ulab_vectorise_module;
 } while(0)
 
 #define MATH_FUN_1(py_name, c_name) \
-    mp_obj_t vectorise_ ## py_name(mp_obj_t x_obj) { \
+    static mp_obj_t vectorise_ ## py_name(mp_obj_t x_obj) { \
         return vectorise_generic_vector(x_obj, MICROPY_FLOAT_C_FUN(c_name)); \
 }
     

docs/manual/source/conf.py

@@ -18,12 +18,11 @@
 # -- Project information -----------------------------------------------------
 
 project = 'micropython-ulab'
-copyright = '2019, Zoltán Vörös'
+copyright = '2019-2020, Zoltán Vörös'
 author = 'Zoltán Vörös'
 
 # The full version, including alpha/beta/rc tags
-release = '0.32'
-
+release = '0.50.0'
 
 # -- General configuration ---------------------------------------------------
 

docs/ulab-change-log.md

@@ -1,3 +1,145 @@
+Mon, 1 Jun 2020
+
+version 0.50.1
+
+    fixes error in numerical max/min
+    
+Mon, 18 May 2020
+
+version 0.50.0
+
+    move interp to the approx sub-module
+    
+Wed, 06 May 2020
+
+version 0.46.0
+
+    add curve_fit to the approx sub-module
+
+version 0.44.0
+
+	add approx sub-module with newton, fmin, and bisect functions
+
+Thu, 30 Apr 2020
+
+version 0.44.0
+
+	add approx sub-module with newton, fmin, and bisect functions
+
+=======
+Tue, 19 May 2020
+
+version 0.46.1
+
+    fixed bad error in binary_op
+
+Wed, 6 May 2020
+
+version 0.46
+
+    added vectorisation of python functions
+
+Sat, 2 May 2020
+
+version 0.45.0
+
+	add equal/not_equal to the compare module
+	
+Tue, 21 Apr 2020
+
+version 0.42.0
+
+	add minimum/maximum/clip functions
+
+Mon, 20 Apr 2020
+
+version 0.41.6
+
+	argument handling improvement in polyfit
+
+Mon, 20 Apr 2020
+
+version 0.41.5
+
+	fix compilation errors due to https://github.com/micropython/micropython/commit/30840ebc9925bb8ef025dbc2d5982b1bfeb75f1b
+
+Sat, 18 Apr 2020
+
+version 0.41.4
+
+	fix compilation error on hardware ports
+
+Tue, 14 Apr 2020
+
+version 0.41.3
+
+	fix indexing error in dot function
+
+Thu, 9 Apr 2020
+
+version 0.41.2
+
+	fix transpose function
+
+Tue, 7 Apr 2020
+
+version 0.41.2
+
+	fix discrepancy in argmin/argmax behaviour
+
+Tue, 7 Apr 2020
+
+version 0.41.1
+
+	fix error in argsort
+
+Sat, 4 Apr 2020
+
+version 0.41.0
+
+	implemented == and != binary operators
+
+Fri, 3 Apr 2020
+
+version 0.40.0
+
+	added trace to linalg
+
+Thu, 2 Apr 2020
+
+version 0.39.0
+
+	added the ** operator, and operand swapping in binary operators
+
+Thu, 2 Apr 2020
+
+version 0.38.1
+
+	added fast option, when initialising from ndarray_properties
+
+Thu, 12 Mar 2020
+
+version 0.38.0
+
+	added initialisation from ndarray, and the around function
+
+Tue, 10 Mar 2020
+
+version 0.37.0
+
+	added Cholesky decomposition to linalg.c
+	
+Thu, 27 Feb 2020
+
+version 0.36.0
+
+	moved zeros, ones, eye and linspace into separate module (they are still bound at the top level)
+	
+Thu, 27 Feb 2020
+
+version 0.35.0
+
+    Move zeros, ones back into top level ulab module
 
 Tue, 18 Feb 2020
 

tests/00smoke.py

@@ -1,2 +1,2 @@
-from ulab import linalg
-print(linalg.eye(3))
+import ulab
+print(ulab.eye(3))

tests/argminmax.py

@@ -0,0 +1,62 @@
+import ulab
+
+# Adapted from https://docs.python.org/3.8/library/itertools.html#itertools.permutations
+def permutations(iterable, r=None):
+    # permutations('ABCD', 2) --> AB AC AD BA BC BD CA CB CD DA DB DC
+    # permutations(range(3)) --> 012 021 102 120 201 210
+    pool = tuple(iterable)
+    n = len(pool)
+    r = n if r is None else r
+    if r > n:
+        return
+    indices = list(range(n))
+    cycles = list(range(n, n-r, -1))
+    yield tuple(pool[i] for i in indices[:r])
+    while n:
+        for i in reversed(range(r)):
+            cycles[i] -= 1
+            if cycles[i] == 0:
+                indices[i:] = indices[i+1:] + indices[i:i+1]
+                cycles[i] = n - i
+            else:
+                j = cycles[i]
+                indices[i], indices[-j] = indices[-j], indices[i]
+                yield tuple(pool[i] for i in indices[:r])
+                break
+        else:
+            return
+
+# Combinations expected to throw
+try:
+    print(ulab.numerical.argmin([]))
+except ValueError:
+    print("ValueError")
+
+try:
+    print(ulab.numerical.argmax([]))
+except ValueError:
+    print("ValueError")
+
+# Combinations expected to succeed
+print(ulab.numerical.argmin([1]))
+print(ulab.numerical.argmax([1]))
+print(ulab.numerical.argmin(ulab.array([1])))
+print(ulab.numerical.argmax(ulab.array([1])))
+
+print()
+print("max tests")
+for p in permutations((100,200,300)):
+    m1 = ulab.numerical.argmax(p)
+    m2 = ulab.numerical.argmax(ulab.array(p))
+    print(p, m1, m2)
+    if m1 != m2 or p[m1] != max(p):
+        print("FAIL", p, m1, m2, max(p))
+
+print()
+print("min tests")
+for p in permutations((100,200,300)):
+    m1 = ulab.numerical.argmin(p)
+    m2 = ulab.numerical.argmin(ulab.array(p))
+    print(p, m1, m2)
+    if m1 != m2 or p[m1] != min(p):
+        print("FAIL", p, m1, m2, min(p))

tests/argminmax.py.exp

@@ -0,0 +1,22 @@
+ValueError
+ValueError
+0
+0
+0
+0
+
+max tests
+(100, 200, 300) 2 2
+(100, 300, 200) 1 1
+(200, 100, 300) 2 2
+(200, 300, 100) 1 1
+(300, 100, 200) 0 0
+(300, 200, 100) 0 0
+
+min tests
+(100, 200, 300) 0 0
+(100, 300, 200) 0 0
+(200, 100, 300) 1 1
+(200, 300, 100) 2 2
+(300, 100, 200) 1 1
+(300, 200, 100) 2 2

tests/cholesky.py

@@ -0,0 +1,17 @@
+import ulab
+from ulab import linalg
+
+a = ulab.array([[1, 2], [2, 5]])
+print(linalg.cholesky(a))
+
+b = a = ulab.array([[25, 15, -5], [15, 18,  0], [-5,  0, 11]])
+print(linalg.cholesky(b))
+
+c = ulab.array([[18, 22,  54,  42], [22, 70,  86,  62], [54, 86, 174, 134], [42, 62, 134, 106]])
+print(linalg.cholesky(c))
+
+# this throw a ValueError exception
+d = ulab.array([[25, 15, -5], [15, 18,  0], [-5,  0, 1]])
+
+
+

tests/cholesky.py.exp

@@ -0,0 +1,9 @@
+array([[1.0, 0.0],
+	 [2.0, 1.0]], dtype=float)
+array([[5.0, 0.0, 0.0],
+	 [3.0, 3.0, 0.0],
+	 [-1.0, 1.0, 3.0]], dtype=float)
+array([[4.242640687119285, 0.0, 0.0, 0.0],
+	 [5.185449728701349, 6.565905201197403, 0.0, 0.0],
+	 [12.72792206135786, 3.046038495400855, 1.649742247909068, 0.0],
+	 [9.899494936611665, 1.624553864213789, 1.849711005231386, 1.392621247645583]], dtype=float)

tests/compare.py

@@ -0,0 +1,14 @@
+import ulab
+from ulab import compare
+
+a = ulab.array([1, 2, 3, 4, 5], dtype=ulab.uint8)
+b = ulab.array([5, 4, 3, 2, 1], dtype=ulab.float)
+print(compare.minimum(a, b))
+print(compare.maximum(a, b))
+print(compare.maximum(1, 5.5))
+
+a = ulab.array(range(9), dtype=ulab.uint8)
+print(compare.clip(a, 3, 7))
+
+b = 3 * ulab.ones(len(a), dtype=ulab.float)
+print(compare.clip(a, b, 7))

tests/compare.py.exp

@@ -0,0 +1,5 @@
+array([1.0, 2.0, 3.0, 2.0, 1.0], dtype=float)
+array([5.0, 4.0, 3.0, 4.0, 5.0], dtype=float)
+5.5
+array([3, 3, 3, 3, 4, 5, 6, 7, 7], dtype=uint8)
+array([3.0, 3.0, 3.0, 3.0, 4.0, 5.0, 6.0, 7.0, 7.0], dtype=float)

tests/constructors.py

@@ -0,0 +1,13 @@
+from ulab import linalg
+import ulab
+print(ulab.ones(3))
+print(ulab.ones((2,3)))
+print(ulab.zeros(3))
+print(ulab.zeros((2,3)))
+print(ulab.eye(3))
+print(ulab.ones(1, dtype=ulab.int8))
+print(ulab.ones(2, dtype=ulab.uint8))
+print(ulab.ones(3, dtype=ulab.int16))
+print(ulab.ones(4, dtype=ulab.uint16))
+print(ulab.ones(5, dtype=ulab.float))
+print(ulab.linspace(0, 1, 9))

tests/constructors.py.exp

@@ -0,0 +1,15 @@
+array([1.0, 1.0, 1.0], dtype=float)
+array([[1.0, 1.0, 1.0],
+	 [1.0, 1.0, 1.0]], dtype=float)
+array([0.0, 0.0, 0.0], dtype=float)
+array([[0.0, 0.0, 0.0],
+	 [0.0, 0.0, 0.0]], dtype=float)
+array([[1.0, 0.0, 0.0],
+	 [0.0, 1.0, 0.0],
+	 [0.0, 0.0, 1.0]], dtype=float)
+array([1], dtype=int8)
+array([1, 1], dtype=uint8)
+array([1, 1, 1], dtype=int16)
+array([1, 1, 1, 1], dtype=uint16)
+array([1.0, 1.0, 1.0, 1.0, 1.0], dtype=float)
+array([0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0], dtype=float)

tests/linalg.py

@@ -0,0 +1,17 @@
+import ulab
+from ulab import linalg
+
+a = ulab.array([[1, 2], [3, 4]])
+print(linalg.inv(a))
+
+b = ulab.array([[1, 2, 3], [4, 5, 6], [7, 8, 7]])
+print(linalg.inv(b))
+
+c = ulab.array([[1, 2, 0, 0], [0, 6, 7, 0], [0, 0, 8, 9], [0, 0, 15, 13]])
+print(linalg.inv(c))
+
+print(linalg.det(a))
+print(linalg.det(b))
+print(linalg.det(c))
+
+

tests/linalg.py.exp

@@ -0,0 +1,12 @@
+array([[-2.0, 1.0],
+	 [1.5, -0.5]], dtype=float)
+array([[-2.166666666666667, 1.666666666666667, -0.5],
+	 [2.333333333333333, -2.333333333333333, 1.0],
+	 [-0.5, 1.0, -0.5]], dtype=float)
+array([[1.0, -0.3333333333333333, -0.9784946236559136, 0.6774193548387095],
+	 [0.0, 0.1666666666666667, 0.489247311827957, -0.3387096774193548],
+	 [0.0, 0.0, -0.4193548387096775, 0.2903225806451613],
+	 [-0.0, -0.0, 0.4838709677419355, -0.2580645161290323]], dtype=float)
+-2.0
+6.0
+-186.0

tests/operators.py

@@ -0,0 +1,20 @@
+import ulab
+a = ulab.ones(3)
+print(a+a)
+print(a-a)
+print(a*a)
+print(a/a)
+print(a+2)
+print(a-2)
+print(a*2)
+print(a/2)
+print(a<1)
+print(a<2)
+print(a<=0)
+print(a<=1)
+print(a>1)
+print(a>2)
+print(a>=0)
+print(a>=1)
+#print(a==0)  # These print just true or false.  Is it right?  is it a micropython limitation?
+#print(a==1)

tests/operators.py.exp

@@ -0,0 +1,16 @@
+array([2.0, 2.0, 2.0], dtype=float)
+array([0.0, 0.0, 0.0], dtype=float)
+array([1.0, 1.0, 1.0], dtype=float)
+array([1.0, 1.0, 1.0], dtype=float)
+array([3.0, 3.0, 3.0], dtype=float)
+array([-1.0, -1.0, -1.0], dtype=float)
+array([2.0, 2.0, 2.0], dtype=float)
+array([0.5, 0.5, 0.5], dtype=float)
+[False, False, False]
+[True, True, True]
+[False, False, False]
+[True, True, True]
+[False, False, False]
+[False, False, False]
+[True, True, True]
+[True, True, True]

tests/poly.py

@@ -0,0 +1,27 @@
+import ulab
+from ulab import poly
+from ulab import vector
+
+# polynom evaluation
+x = ulab.linspace(0, 10, num=9)
+p = [1, 2, 3]
+y = poly.polyval(p, x)
+print(y)
+
+# linear fit
+x = ulab.linspace(-5, 5, num=11)
+y = x + vector.sin(x)
+p = poly.polyfit(x, y, 1)
+print(p)
+
+# quadratic fit
+x = ulab.linspace(-5, 5, num=11)
+y = x*x + vector.sin(x)*3.0
+p = poly.polyfit(x, y, 2)
+print(p)
+
+# cubic fit
+x = ulab.linspace(-5, 5, num=11)
+y = x*x*x + vector.sin(x)*10.0
+p = poly.polyfit(x, y, 3)
+print(p)

tests/poly.py.exp

@@ -0,0 +1,4 @@
+array([3.0, 7.0625, 14.25, 24.5625, 38.0, 54.56250000000001, 74.25000000000001, 97.06250000000001, 123.0], dtype=float)
+array([0.9138471728743898, -8.074349270001139e-17], dtype=float)
+array([1.0, -0.2584584813768293, 3.552713678800501e-15], dtype=float)
+array([0.766798803225857, 0.0, 3.289453031323674, 0.0], dtype=float)

tests/slicing.py

@@ -0,0 +1,25 @@
+try:
+    import ulab as np
+except:
+    import numpy as np
+
+for num in range(1,4):
+    for start in range(-num, num+1):
+        for end in range(-num, num+1):
+            for stride in (-3, -2, -1, 1, 2, 3):
+                l = list(range(num))
+                a = np.array(l, dtype=np.int8)
+                sl = l[start:end:stride]
+                ll = len(sl)
+                try:
+                    sa = list(a[start:end:stride])
+                    la = len(sa)
+                except IndexError as e:
+                    sa = str(e)
+                    la = -1
+                print("%2d [% d:% d:% d]     %-24r %-24r%s" % (
+                    num, start, end, stride, sl, sa, " ***" if sa != sl else ""))
+
+                a[start:end:stride] = np.ones(len(sl)) * -1
+                print("%2d [% d:% d:% d]     %r" % (
+                    num, start, end, stride, list(a)))

tests/slicing.py.exp

@@ -0,0 +1,996 @@
+ 1 [-1:-1:-3]     []                       []                      
+ 1 [-1:-1:-3]     [0]
+ 1 [-1:-1:-2]     []                       []                      
+ 1 [-1:-1:-2]     [0]
+ 1 [-1:-1:-1]     []                       []                      
+ 1 [-1:-1:-1]     [0]
+ 1 [-1:-1: 1]     []                       []                      
+ 1 [-1:-1: 1]     [0]
+ 1 [-1:-1: 2]     []                       []                      
+ 1 [-1:-1: 2]     [0]
+ 1 [-1:-1: 3]     []                       []                      
+ 1 [-1:-1: 3]     [0]
+ 1 [-1: 0:-3]     []                       []                      
+ 1 [-1: 0:-3]     [0]
+ 1 [-1: 0:-2]     []                       []                      
+ 1 [-1: 0:-2]     [0]
+ 1 [-1: 0:-1]     []                       []                      
+ 1 [-1: 0:-1]     [0]
+ 1 [-1: 0: 1]     []                       []                      
+ 1 [-1: 0: 1]     [0]
+ 1 [-1: 0: 2]     []                       []                      
+ 1 [-1: 0: 2]     [0]
+ 1 [-1: 0: 3]     []                       []                      
+ 1 [-1: 0: 3]     [0]
+ 1 [-1: 1:-3]     []                       []                      
+ 1 [-1: 1:-3]     [0]
+ 1 [-1: 1:-2]     []                       []                      
+ 1 [-1: 1:-2]     [0]
+ 1 [-1: 1:-1]     []                       []                      
+ 1 [-1: 1:-1]     [0]
+ 1 [-1: 1: 1]     [0]                      [0]                     
+ 1 [-1: 1: 1]     [-1]
+ 1 [-1: 1: 2]     [0]                      [0]                     
+ 1 [-1: 1: 2]     [-1]
+ 1 [-1: 1: 3]     [0]                      [0]                     
+ 1 [-1: 1: 3]     [-1]
+ 1 [ 0:-1:-3]     []                       []                      
+ 1 [ 0:-1:-3]     [0]
+ 1 [ 0:-1:-2]     []                       []                      
+ 1 [ 0:-1:-2]     [0]
+ 1 [ 0:-1:-1]     []                       []                      
+ 1 [ 0:-1:-1]     [0]
+ 1 [ 0:-1: 1]     []                       []                      
+ 1 [ 0:-1: 1]     [0]
+ 1 [ 0:-1: 2]     []                       []                      
+ 1 [ 0:-1: 2]     [0]
+ 1 [ 0:-1: 3]     []                       []                      
+ 1 [ 0:-1: 3]     [0]
+ 1 [ 0: 0:-3]     []                       []                      
+ 1 [ 0: 0:-3]     [0]
+ 1 [ 0: 0:-2]     []                       []                      
+ 1 [ 0: 0:-2]     [0]
+ 1 [ 0: 0:-1]     []                       []                      
+ 1 [ 0: 0:-1]     [0]
+ 1 [ 0: 0: 1]     []                       []                      
+ 1 [ 0: 0: 1]     [0]
+ 1 [ 0: 0: 2]     []                       []                      
+ 1 [ 0: 0: 2]     [0]
+ 1 [ 0: 0: 3]     []                       []                      
+ 1 [ 0: 0: 3]     [0]
+ 1 [ 0: 1:-3]     []                       []                      
+ 1 [ 0: 1:-3]     [0]
+ 1 [ 0: 1:-2]     []                       []                      
+ 1 [ 0: 1:-2]     [0]
+ 1 [ 0: 1:-1]     []                       []                      
+ 1 [ 0: 1:-1]     [0]
+ 1 [ 0: 1: 1]     [0]                      [0]                     
+ 1 [ 0: 1: 1]     [-1]
+ 1 [ 0: 1: 2]     [0]                      [0]                     
+ 1 [ 0: 1: 2]     [-1]
+ 1 [ 0: 1: 3]     [0]                      [0]                     
+ 1 [ 0: 1: 3]     [-1]
+ 1 [ 1:-1:-3]     []                       []                      
+ 1 [ 1:-1:-3]     [0]
+ 1 [ 1:-1:-2]     []                       []                      
+ 1 [ 1:-1:-2]     [0]
+ 1 [ 1:-1:-1]     []                       []                      
+ 1 [ 1:-1:-1]     [0]
+ 1 [ 1:-1: 1]     []                       []                      
+ 1 [ 1:-1: 1]     [0]
+ 1 [ 1:-1: 2]     []                       []                      
+ 1 [ 1:-1: 2]     [0]
+ 1 [ 1:-1: 3]     []                       []                      
+ 1 [ 1:-1: 3]     [0]
+ 1 [ 1: 0:-3]     []                       []                      
+ 1 [ 1: 0:-3]     [0]
+ 1 [ 1: 0:-2]     []                       []                      
+ 1 [ 1: 0:-2]     [0]
+ 1 [ 1: 0:-1]     []                       []                      
+ 1 [ 1: 0:-1]     [0]
+ 1 [ 1: 0: 1]     []                       []                      
+ 1 [ 1: 0: 1]     [0]
+ 1 [ 1: 0: 2]     []                       []                      
+ 1 [ 1: 0: 2]     [0]
+ 1 [ 1: 0: 3]     []                       []                      
+ 1 [ 1: 0: 3]     [0]
+ 1 [ 1: 1:-3]     []                       []                      
+ 1 [ 1: 1:-3]     [0]
+ 1 [ 1: 1:-2]     []                       []                      
+ 1 [ 1: 1:-2]     [0]
+ 1 [ 1: 1:-1]     []                       []                      
+ 1 [ 1: 1:-1]     [0]
+ 1 [ 1: 1: 1]     []                       []                      
+ 1 [ 1: 1: 1]     [0]
+ 1 [ 1: 1: 2]     []                       []                      
+ 1 [ 1: 1: 2]     [0]
+ 1 [ 1: 1: 3]     []                       []                      
+ 1 [ 1: 1: 3]     [0]
+ 2 [-2:-2:-3]     []                       []                      
+ 2 [-2:-2:-3]     [0, 1]
+ 2 [-2:-2:-2]     []                       []                      
+ 2 [-2:-2:-2]     [0, 1]
+ 2 [-2:-2:-1]     []                       []                      
+ 2 [-2:-2:-1]     [0, 1]
+ 2 [-2:-2: 1]     []                       []                      
+ 2 [-2:-2: 1]     [0, 1]
+ 2 [-2:-2: 2]     []                       []                      
+ 2 [-2:-2: 2]     [0, 1]
+ 2 [-2:-2: 3]     []                       []                      
+ 2 [-2:-2: 3]     [0, 1]
+ 2 [-2:-1:-3]     []                       []                      
+ 2 [-2:-1:-3]     [0, 1]
+ 2 [-2:-1:-2]     []                       []                      
+ 2 [-2:-1:-2]     [0, 1]
+ 2 [-2:-1:-1]     []                       []                      
+ 2 [-2:-1:-1]     [0, 1]
+ 2 [-2:-1: 1]     [0]                      [0]                     
+ 2 [-2:-1: 1]     [-1, 1]
+ 2 [-2:-1: 2]     [0]                      [0]                     
+ 2 [-2:-1: 2]     [-1, 1]
+ 2 [-2:-1: 3]     [0]                      [0]                     
+ 2 [-2:-1: 3]     [-1, 1]
+ 2 [-2: 0:-3]     []                       []                      
+ 2 [-2: 0:-3]     [0, 1]
+ 2 [-2: 0:-2]     []                       []                      
+ 2 [-2: 0:-2]     [0, 1]
+ 2 [-2: 0:-1]     []                       []                      
+ 2 [-2: 0:-1]     [0, 1]
+ 2 [-2: 0: 1]     []                       []                      
+ 2 [-2: 0: 1]     [0, 1]
+ 2 [-2: 0: 2]     []                       []                      
+ 2 [-2: 0: 2]     [0, 1]
+ 2 [-2: 0: 3]     []                       []                      
+ 2 [-2: 0: 3]     [0, 1]
+ 2 [-2: 1:-3]     []                       []                      
+ 2 [-2: 1:-3]     [0, 1]
+ 2 [-2: 1:-2]     []                       []                      
+ 2 [-2: 1:-2]     [0, 1]
+ 2 [-2: 1:-1]     []                       []                      
+ 2 [-2: 1:-1]     [0, 1]
+ 2 [-2: 1: 1]     [0]                      [0]                     
+ 2 [-2: 1: 1]     [-1, 1]
+ 2 [-2: 1: 2]     [0]                      [0]                     
+ 2 [-2: 1: 2]     [-1, 1]
+ 2 [-2: 1: 3]     [0]                      [0]                     
+ 2 [-2: 1: 3]     [-1, 1]
+ 2 [-2: 2:-3]     []                       []                      
+ 2 [-2: 2:-3]     [0, 1]
+ 2 [-2: 2:-2]     []                       []                      
+ 2 [-2: 2:-2]     [0, 1]
+ 2 [-2: 2:-1]     []                       []                      
+ 2 [-2: 2:-1]     [0, 1]
+ 2 [-2: 2: 1]     [0, 1]                   [0, 1]                  
+ 2 [-2: 2: 1]     [-1, -1]
+ 2 [-2: 2: 2]     [0]                      [0]                     
+ 2 [-2: 2: 2]     [-1, 1]
+ 2 [-2: 2: 3]     [0]                      [0]                     
+ 2 [-2: 2: 3]     [-1, 1]
+ 2 [-1:-2:-3]     [1]                      [1]                     
+ 2 [-1:-2:-3]     [0, -1]
+ 2 [-1:-2:-2]     [1]                      [1]                     
+ 2 [-1:-2:-2]     [0, -1]
+ 2 [-1:-2:-1]     [1]                      [1]                     
+ 2 [-1:-2:-1]     [0, -1]
+ 2 [-1:-2: 1]     []                       []                      
+ 2 [-1:-2: 1]     [0, 1]
+ 2 [-1:-2: 2]     []                       []                      
+ 2 [-1:-2: 2]     [0, 1]
+ 2 [-1:-2: 3]     []                       []                      
+ 2 [-1:-2: 3]     [0, 1]
+ 2 [-1:-1:-3]     []                       []                      
+ 2 [-1:-1:-3]     [0, 1]
+ 2 [-1:-1:-2]     []                       []                      
+ 2 [-1:-1:-2]     [0, 1]
+ 2 [-1:-1:-1]     []                       []                      
+ 2 [-1:-1:-1]     [0, 1]
+ 2 [-1:-1: 1]     []                       []                      
+ 2 [-1:-1: 1]     [0, 1]
+ 2 [-1:-1: 2]     []                       []                      
+ 2 [-1:-1: 2]     [0, 1]
+ 2 [-1:-1: 3]     []                       []                      
+ 2 [-1:-1: 3]     [0, 1]
+ 2 [-1: 0:-3]     [1]                      [1]                     
+ 2 [-1: 0:-3]     [0, -1]
+ 2 [-1: 0:-2]     [1]                      [1]                     
+ 2 [-1: 0:-2]     [0, -1]
+ 2 [-1: 0:-1]     [1]                      [1]                     
+ 2 [-1: 0:-1]     [0, -1]
+ 2 [-1: 0: 1]     []                       []                      
+ 2 [-1: 0: 1]     [0, 1]
+ 2 [-1: 0: 2]     []                       []                      
+ 2 [-1: 0: 2]     [0, 1]
+ 2 [-1: 0: 3]     []                       []                      
+ 2 [-1: 0: 3]     [0, 1]
+ 2 [-1: 1:-3]     []                       []                      
+ 2 [-1: 1:-3]     [0, 1]
+ 2 [-1: 1:-2]     []                       []                      
+ 2 [-1: 1:-2]     [0, 1]
+ 2 [-1: 1:-1]     []                       []                      
+ 2 [-1: 1:-1]     [0, 1]
+ 2 [-1: 1: 1]     []                       []                      
+ 2 [-1: 1: 1]     [0, 1]
+ 2 [-1: 1: 2]     []                       []                      
+ 2 [-1: 1: 2]     [0, 1]
+ 2 [-1: 1: 3]     []                       []                      
+ 2 [-1: 1: 3]     [0, 1]
+ 2 [-1: 2:-3]     []                       []                      
+ 2 [-1: 2:-3]     [0, 1]
+ 2 [-1: 2:-2]     []                       []                      
+ 2 [-1: 2:-2]     [0, 1]
+ 2 [-1: 2:-1]     []                       []                      
+ 2 [-1: 2:-1]     [0, 1]
+ 2 [-1: 2: 1]     [1]                      [1]                     
+ 2 [-1: 2: 1]     [0, -1]
+ 2 [-1: 2: 2]     [1]                      [1]                     
+ 2 [-1: 2: 2]     [0, -1]
+ 2 [-1: 2: 3]     [1]                      [1]                     
+ 2 [-1: 2: 3]     [0, -1]
+ 2 [ 0:-2:-3]     []                       []                      
+ 2 [ 0:-2:-3]     [0, 1]
+ 2 [ 0:-2:-2]     []                       []                      
+ 2 [ 0:-2:-2]     [0, 1]
+ 2 [ 0:-2:-1]     []                       []                      
+ 2 [ 0:-2:-1]     [0, 1]
+ 2 [ 0:-2: 1]     []                       []                      
+ 2 [ 0:-2: 1]     [0, 1]
+ 2 [ 0:-2: 2]     []                       []                      
+ 2 [ 0:-2: 2]     [0, 1]
+ 2 [ 0:-2: 3]     []                       []                      
+ 2 [ 0:-2: 3]     [0, 1]
+ 2 [ 0:-1:-3]     []                       []                      
+ 2 [ 0:-1:-3]     [0, 1]
+ 2 [ 0:-1:-2]     []                       []                      
+ 2 [ 0:-1:-2]     [0, 1]
+ 2 [ 0:-1:-1]     []                       []                      
+ 2 [ 0:-1:-1]     [0, 1]
+ 2 [ 0:-1: 1]     [0]                      [0]                     
+ 2 [ 0:-1: 1]     [-1, 1]
+ 2 [ 0:-1: 2]     [0]                      [0]                     
+ 2 [ 0:-1: 2]     [-1, 1]
+ 2 [ 0:-1: 3]     [0]                      [0]                     
+ 2 [ 0:-1: 3]     [-1, 1]
+ 2 [ 0: 0:-3]     []                       []                      
+ 2 [ 0: 0:-3]     [0, 1]
+ 2 [ 0: 0:-2]     []                       []                      
+ 2 [ 0: 0:-2]     [0, 1]
+ 2 [ 0: 0:-1]     []                       []                      
+ 2 [ 0: 0:-1]     [0, 1]
+ 2 [ 0: 0: 1]     []                       []                      
+ 2 [ 0: 0: 1]     [0, 1]
+ 2 [ 0: 0: 2]     []                       []                      
+ 2 [ 0: 0: 2]     [0, 1]
+ 2 [ 0: 0: 3]     []                       []                      
+ 2 [ 0: 0: 3]     [0, 1]
+ 2 [ 0: 1:-3]     []                       []                      
+ 2 [ 0: 1:-3]     [0, 1]
+ 2 [ 0: 1:-2]     []                       []                      
+ 2 [ 0: 1:-2]     [0, 1]
+ 2 [ 0: 1:-1]     []                       []                      
+ 2 [ 0: 1:-1]     [0, 1]
+ 2 [ 0: 1: 1]     [0]                      [0]                     
+ 2 [ 0: 1: 1]     [-1, 1]
+ 2 [ 0: 1: 2]     [0]                      [0]                     
+ 2 [ 0: 1: 2]     [-1, 1]
+ 2 [ 0: 1: 3]     [0]                      [0]                     
+ 2 [ 0: 1: 3]     [-1, 1]
+ 2 [ 0: 2:-3]     []                       []                      
+ 2 [ 0: 2:-3]     [0, 1]
+ 2 [ 0: 2:-2]     []                       []                      
+ 2 [ 0: 2:-2]     [0, 1]
+ 2 [ 0: 2:-1]     []                       []                      
+ 2 [ 0: 2:-1]     [0, 1]
+ 2 [ 0: 2: 1]     [0, 1]                   [0, 1]                  
+ 2 [ 0: 2: 1]     [-1, -1]
+ 2 [ 0: 2: 2]     [0]                      [0]                     
+ 2 [ 0: 2: 2]     [-1, 1]
+ 2 [ 0: 2: 3]     [0]                      [0]                     
+ 2 [ 0: 2: 3]     [-1, 1]
+ 2 [ 1:-2:-3]     [1]                      [1]                     
+ 2 [ 1:-2:-3]     [0, -1]
+ 2 [ 1:-2:-2]     [1]                      [1]                     
+ 2 [ 1:-2:-2]     [0, -1]
+ 2 [ 1:-2:-1]     [1]                      [1]                     
+ 2 [ 1:-2:-1]     [0, -1]
+ 2 [ 1:-2: 1]     []                       []                      
+ 2 [ 1:-2: 1]     [0, 1]
+ 2 [ 1:-2: 2]     []                       []                      
+ 2 [ 1:-2: 2]     [0, 1]
+ 2 [ 1:-2: 3]     []                       []                      
+ 2 [ 1:-2: 3]     [0, 1]
+ 2 [ 1:-1:-3]     []                       []                      
+ 2 [ 1:-1:-3]     [0, 1]
+ 2 [ 1:-1:-2]     []                       []                      
+ 2 [ 1:-1:-2]     [0, 1]
+ 2 [ 1:-1:-1]     []                       []                      
+ 2 [ 1:-1:-1]     [0, 1]
+ 2 [ 1:-1: 1]     []                       []                      
+ 2 [ 1:-1: 1]     [0, 1]
+ 2 [ 1:-1: 2]     []                       []                      
+ 2 [ 1:-1: 2]     [0, 1]
+ 2 [ 1:-1: 3]     []                       []                      
+ 2 [ 1:-1: 3]     [0, 1]
+ 2 [ 1: 0:-3]     [1]                      [1]                     
+ 2 [ 1: 0:-3]     [0, -1]
+ 2 [ 1: 0:-2]     [1]                      [1]                     
+ 2 [ 1: 0:-2]     [0, -1]
+ 2 [ 1: 0:-1]     [1]                      [1]                     
+ 2 [ 1: 0:-1]     [0, -1]
+ 2 [ 1: 0: 1]     []                       []                      
+ 2 [ 1: 0: 1]     [0, 1]
+ 2 [ 1: 0: 2]     []                       []                      
+ 2 [ 1: 0: 2]     [0, 1]
+ 2 [ 1: 0: 3]     []                       []                      
+ 2 [ 1: 0: 3]     [0, 1]
+ 2 [ 1: 1:-3]     []                       []                      
+ 2 [ 1: 1:-3]     [0, 1]
+ 2 [ 1: 1:-2]     []                       []                      
+ 2 [ 1: 1:-2]     [0, 1]
+ 2 [ 1: 1:-1]     []                       []                      
+ 2 [ 1: 1:-1]     [0, 1]
+ 2 [ 1: 1: 1]     []                       []                      
+ 2 [ 1: 1: 1]     [0, 1]
+ 2 [ 1: 1: 2]     []                       []                      
+ 2 [ 1: 1: 2]     [0, 1]
+ 2 [ 1: 1: 3]     []                       []                      
+ 2 [ 1: 1: 3]     [0, 1]
+ 2 [ 1: 2:-3]     []                       []                      
+ 2 [ 1: 2:-3]     [0, 1]
+ 2 [ 1: 2:-2]     []                       []                      
+ 2 [ 1: 2:-2]     [0, 1]
+ 2 [ 1: 2:-1]     []                       []                      
+ 2 [ 1: 2:-1]     [0, 1]
+ 2 [ 1: 2: 1]     [1]                      [1]                     
+ 2 [ 1: 2: 1]     [0, -1]
+ 2 [ 1: 2: 2]     [1]                      [1]                     
+ 2 [ 1: 2: 2]     [0, -1]
+ 2 [ 1: 2: 3]     [1]                      [1]                     
+ 2 [ 1: 2: 3]     [0, -1]
+ 2 [ 2:-2:-3]     [1]                      [1]                     
+ 2 [ 2:-2:-3]     [0, -1]
+ 2 [ 2:-2:-2]     [1]                      [1]                     
+ 2 [ 2:-2:-2]     [0, -1]
+ 2 [ 2:-2:-1]     [1]                      [1]                     
+ 2 [ 2:-2:-1]     [0, -1]
+ 2 [ 2:-2: 1]     []                       []                      
+ 2 [ 2:-2: 1]     [0, 1]
+ 2 [ 2:-2: 2]     []                       []                      
+ 2 [ 2:-2: 2]     [0, 1]
+ 2 [ 2:-2: 3]     []                       []                      
+ 2 [ 2:-2: 3]     [0, 1]
+ 2 [ 2:-1:-3]     []                       []                      
+ 2 [ 2:-1:-3]     [0, 1]
+ 2 [ 2:-1:-2]     []                       []                      
+ 2 [ 2:-1:-2]     [0, 1]
+ 2 [ 2:-1:-1]     []                       []                      
+ 2 [ 2:-1:-1]     [0, 1]
+ 2 [ 2:-1: 1]     []                       []                      
+ 2 [ 2:-1: 1]     [0, 1]
+ 2 [ 2:-1: 2]     []                       []                      
+ 2 [ 2:-1: 2]     [0, 1]
+ 2 [ 2:-1: 3]     []                       []                      
+ 2 [ 2:-1: 3]     [0, 1]
+ 2 [ 2: 0:-3]     [1]                      [1]                     
+ 2 [ 2: 0:-3]     [0, -1]
+ 2 [ 2: 0:-2]     [1]                      [1]                     
+ 2 [ 2: 0:-2]     [0, -1]
+ 2 [ 2: 0:-1]     [1]                      [1]                     
+ 2 [ 2: 0:-1]     [0, -1]
+ 2 [ 2: 0: 1]     []                       []                      
+ 2 [ 2: 0: 1]     [0, 1]
+ 2 [ 2: 0: 2]     []                       []                      
+ 2 [ 2: 0: 2]     [0, 1]
+ 2 [ 2: 0: 3]     []                       []                      
+ 2 [ 2: 0: 3]     [0, 1]
+ 2 [ 2: 1:-3]     []                       []                      
+ 2 [ 2: 1:-3]     [0, 1]
+ 2 [ 2: 1:-2]     []                       []                      
+ 2 [ 2: 1:-2]     [0, 1]
+ 2 [ 2: 1:-1]     []                       []                      
+ 2 [ 2: 1:-1]     [0, 1]
+ 2 [ 2: 1: 1]     []                       []                      
+ 2 [ 2: 1: 1]     [0, 1]
+ 2 [ 2: 1: 2]     []                       []                      
+ 2 [ 2: 1: 2]     [0, 1]
+ 2 [ 2: 1: 3]     []                       []                      
+ 2 [ 2: 1: 3]     [0, 1]
+ 2 [ 2: 2:-3]     []                       []                      
+ 2 [ 2: 2:-3]     [0, 1]
+ 2 [ 2: 2:-2]     []                       []                      
+ 2 [ 2: 2:-2]     [0, 1]
+ 2 [ 2: 2:-1]     []                       []                      
+ 2 [ 2: 2:-1]     [0, 1]
+ 2 [ 2: 2: 1]     []                       []                      
+ 2 [ 2: 2: 1]     [0, 1]
+ 2 [ 2: 2: 2]     []                       []                      
+ 2 [ 2: 2: 2]     [0, 1]
+ 2 [ 2: 2: 3]     []                       []                      
+ 2 [ 2: 2: 3]     [0, 1]
+ 3 [-3:-3:-3]     []                       []                      
+ 3 [-3:-3:-3]     [0, 1, 2]
+ 3 [-3:-3:-2]     []                       []                      
+ 3 [-3:-3:-2]     [0, 1, 2]
+ 3 [-3:-3:-1]     []                       []                      
+ 3 [-3:-3:-1]     [0, 1, 2]
+ 3 [-3:-3: 1]     []                       []                      
+ 3 [-3:-3: 1]     [0, 1, 2]
+ 3 [-3:-3: 2]     []                       []                      
+ 3 [-3:-3: 2]     [0, 1, 2]
+ 3 [-3:-3: 3]     []                       []                      
+ 3 [-3:-3: 3]     [0, 1, 2]
+ 3 [-3:-2:-3]     []                       []                      
+ 3 [-3:-2:-3]     [0, 1, 2]
+ 3 [-3:-2:-2]     []                       []                      
+ 3 [-3:-2:-2]     [0, 1, 2]
+ 3 [-3:-2:-1]     []                       []                      
+ 3 [-3:-2:-1]     [0, 1, 2]
+ 3 [-3:-2: 1]     [0]                      [0]                     
+ 3 [-3:-2: 1]     [-1, 1, 2]
+ 3 [-3:-2: 2]     [0]                      [0]                     
+ 3 [-3:-2: 2]     [-1, 1, 2]
+ 3 [-3:-2: 3]     [0]                      [0]                     
+ 3 [-3:-2: 3]     [-1, 1, 2]
+ 3 [-3:-1:-3]     []                       []                      
+ 3 [-3:-1:-3]     [0, 1, 2]
+ 3 [-3:-1:-2]     []                       []                      
+ 3 [-3:-1:-2]     [0, 1, 2]
+ 3 [-3:-1:-1]     []                       []                      
+ 3 [-3:-1:-1]     [0, 1, 2]
+ 3 [-3:-1: 1]     [0, 1]                   [0, 1]                  
+ 3 [-3:-1: 1]     [-1, -1, 2]
+ 3 [-3:-1: 2]     [0]                      [0]                     
+ 3 [-3:-1: 2]     [-1, 1, 2]
+ 3 [-3:-1: 3]     [0]                      [0]                     
+ 3 [-3:-1: 3]     [-1, 1, 2]
+ 3 [-3: 0:-3]     []                       []                      
+ 3 [-3: 0:-3]     [0, 1, 2]
+ 3 [-3: 0:-2]     []                       []                      
+ 3 [-3: 0:-2]     [0, 1, 2]
+ 3 [-3: 0:-1]     []                       []                      
+ 3 [-3: 0:-1]     [0, 1, 2]
+ 3 [-3: 0: 1]     []                       []                      
+ 3 [-3: 0: 1]     [0, 1, 2]
+ 3 [-3: 0: 2]     []                       []                      
+ 3 [-3: 0: 2]     [0, 1, 2]
+ 3 [-3: 0: 3]     []                       []                      
+ 3 [-3: 0: 3]     [0, 1, 2]
+ 3 [-3: 1:-3]     []                       []                      
+ 3 [-3: 1:-3]     [0, 1, 2]
+ 3 [-3: 1:-2]     []                       []                      
+ 3 [-3: 1:-2]     [0, 1, 2]
+ 3 [-3: 1:-1]     []                       []                      
+ 3 [-3: 1:-1]     [0, 1, 2]
+ 3 [-3: 1: 1]     [0]                      [0]                     
+ 3 [-3: 1: 1]     [-1, 1, 2]
+ 3 [-3: 1: 2]     [0]                      [0]                     
+ 3 [-3: 1: 2]     [-1, 1, 2]
+ 3 [-3: 1: 3]     [0]                      [0]                     
+ 3 [-3: 1: 3]     [-1, 1, 2]
+ 3 [-3: 2:-3]     []                       []                      
+ 3 [-3: 2:-3]     [0, 1, 2]
+ 3 [-3: 2:-2]     []                       []                      
+ 3 [-3: 2:-2]     [0, 1, 2]
+ 3 [-3: 2:-1]     []                       []                      
+ 3 [-3: 2:-1]     [0, 1, 2]
+ 3 [-3: 2: 1]     [0, 1]                   [0, 1]                  
+ 3 [-3: 2: 1]     [-1, -1, 2]
+ 3 [-3: 2: 2]     [0]                      [0]                     
+ 3 [-3: 2: 2]     [-1, 1, 2]
+ 3 [-3: 2: 3]     [0]                      [0]                     
+ 3 [-3: 2: 3]     [-1, 1, 2]
+ 3 [-3: 3:-3]     []                       []                      
+ 3 [-3: 3:-3]     [0, 1, 2]
+ 3 [-3: 3:-2]     []                       []                      
+ 3 [-3: 3:-2]     [0, 1, 2]
+ 3 [-3: 3:-1]     []                       []                      
+ 3 [-3: 3:-1]     [0, 1, 2]
+ 3 [-3: 3: 1]     [0, 1, 2]                [0, 1, 2]               
+ 3 [-3: 3: 1]     [-1, -1, -1]
+ 3 [-3: 3: 2]     [0, 2]                   [0, 2]                  
+ 3 [-3: 3: 2]     [-1, 1, -1]
+ 3 [-3: 3: 3]     [0]                      [0]                     
+ 3 [-3: 3: 3]     [-1, 1, 2]
+ 3 [-2:-3:-3]     [1]                      [1]                     
+ 3 [-2:-3:-3]     [0, -1, 2]
+ 3 [-2:-3:-2]     [1]                      [1]                     
+ 3 [-2:-3:-2]     [0, -1, 2]
+ 3 [-2:-3:-1]     [1]                      [1]                     
+ 3 [-2:-3:-1]     [0, -1, 2]
+ 3 [-2:-3: 1]     []                       []                      
+ 3 [-2:-3: 1]     [0, 1, 2]
+ 3 [-2:-3: 2]     []                       []                      
+ 3 [-2:-3: 2]     [0, 1, 2]
+ 3 [-2:-3: 3]     []                       []                      
+ 3 [-2:-3: 3]     [0, 1, 2]
+ 3 [-2:-2:-3]     []                       []                      
+ 3 [-2:-2:-3]     [0, 1, 2]
+ 3 [-2:-2:-2]     []                       []                      
+ 3 [-2:-2:-2]     [0, 1, 2]
+ 3 [-2:-2:-1]     []                       []                      
+ 3 [-2:-2:-1]     [0, 1, 2]
+ 3 [-2:-2: 1]     []                       []                      
+ 3 [-2:-2: 1]     [0, 1, 2]
+ 3 [-2:-2: 2]     []                       []                      
+ 3 [-2:-2: 2]     [0, 1, 2]
+ 3 [-2:-2: 3]     []                       []                      
+ 3 [-2:-2: 3]     [0, 1, 2]
+ 3 [-2:-1:-3]     []                       []                      
+ 3 [-2:-1:-3]     [0, 1, 2]
+ 3 [-2:-1:-2]     []                       []                      
+ 3 [-2:-1:-2]     [0, 1, 2]
+ 3 [-2:-1:-1]     []                       []                      
+ 3 [-2:-1:-1]     [0, 1, 2]
+ 3 [-2:-1: 1]     [1]                      [1]                     
+ 3 [-2:-1: 1]     [0, -1, 2]
+ 3 [-2:-1: 2]     [1]                      [1]                     
+ 3 [-2:-1: 2]     [0, -1, 2]
+ 3 [-2:-1: 3]     [1]                      [1]                     
+ 3 [-2:-1: 3]     [0, -1, 2]
+ 3 [-2: 0:-3]     [1]                      [1]                     
+ 3 [-2: 0:-3]     [0, -1, 2]
+ 3 [-2: 0:-2]     [1]                      [1]                     
+ 3 [-2: 0:-2]     [0, -1, 2]
+ 3 [-2: 0:-1]     [1]                      [1]                     
+ 3 [-2: 0:-1]     [0, -1, 2]
+ 3 [-2: 0: 1]     []                       []                      
+ 3 [-2: 0: 1]     [0, 1, 2]
+ 3 [-2: 0: 2]     []                       []                      
+ 3 [-2: 0: 2]     [0, 1, 2]
+ 3 [-2: 0: 3]     []                       []                      
+ 3 [-2: 0: 3]     [0, 1, 2]
+ 3 [-2: 1:-3]     []                       []                      
+ 3 [-2: 1:-3]     [0, 1, 2]
+ 3 [-2: 1:-2]     []                       []                      
+ 3 [-2: 1:-2]     [0, 1, 2]
+ 3 [-2: 1:-1]     []                       []                      
+ 3 [-2: 1:-1]     [0, 1, 2]
+ 3 [-2: 1: 1]     []                       []                      
+ 3 [-2: 1: 1]     [0, 1, 2]
+ 3 [-2: 1: 2]     []                       []                      
+ 3 [-2: 1: 2]     [0, 1, 2]
+ 3 [-2: 1: 3]     []                       []                      
+ 3 [-2: 1: 3]     [0, 1, 2]
+ 3 [-2: 2:-3]     []                       []                      
+ 3 [-2: 2:-3]     [0, 1, 2]
+ 3 [-2: 2:-2]     []                       []                      
+ 3 [-2: 2:-2]     [0, 1, 2]
+ 3 [-2: 2:-1]     []                       []                      
+ 3 [-2: 2:-1]     [0, 1, 2]
+ 3 [-2: 2: 1]     [1]                      [1]                     
+ 3 [-2: 2: 1]     [0, -1, 2]
+ 3 [-2: 2: 2]     [1]                      [1]                     
+ 3 [-2: 2: 2]     [0, -1, 2]
+ 3 [-2: 2: 3]     [1]                      [1]                     
+ 3 [-2: 2: 3]     [0, -1, 2]
+ 3 [-2: 3:-3]     []                       []                      
+ 3 [-2: 3:-3]     [0, 1, 2]
+ 3 [-2: 3:-2]     []                       []                      
+ 3 [-2: 3:-2]     [0, 1, 2]
+ 3 [-2: 3:-1]     []                       []                      
+ 3 [-2: 3:-1]     [0, 1, 2]
+ 3 [-2: 3: 1]     [1, 2]                   [1, 2]                  
+ 3 [-2: 3: 1]     [0, -1, -1]
+ 3 [-2: 3: 2]     [1]                      [1]                     
+ 3 [-2: 3: 2]     [0, -1, 2]
+ 3 [-2: 3: 3]     [1]                      [1]                     
+ 3 [-2: 3: 3]     [0, -1, 2]
+ 3 [-1:-3:-3]     [2]                      [2]                     
+ 3 [-1:-3:-3]     [0, 1, -1]
+ 3 [-1:-3:-2]     [2]                      [2]                     
+ 3 [-1:-3:-2]     [0, 1, -1]
+ 3 [-1:-3:-1]     [2, 1]                   [2, 1]                  
+ 3 [-1:-3:-1]     [0, -1, -1]
+ 3 [-1:-3: 1]     []                       []                      
+ 3 [-1:-3: 1]     [0, 1, 2]
+ 3 [-1:-3: 2]     []                       []                      
+ 3 [-1:-3: 2]     [0, 1, 2]
+ 3 [-1:-3: 3]     []                       []                      
+ 3 [-1:-3: 3]     [0, 1, 2]
+ 3 [-1:-2:-3]     [2]                      [2]                     
+ 3 [-1:-2:-3]     [0, 1, -1]
+ 3 [-1:-2:-2]     [2]                      [2]                     
+ 3 [-1:-2:-2]     [0, 1, -1]
+ 3 [-1:-2:-1]     [2]                      [2]                     
+ 3 [-1:-2:-1]     [0, 1, -1]
+ 3 [-1:-2: 1]     []                       []                      
+ 3 [-1:-2: 1]     [0, 1, 2]
+ 3 [-1:-2: 2]     []                       []                      
+ 3 [-1:-2: 2]     [0, 1, 2]
+ 3 [-1:-2: 3]     []                       []                      
+ 3 [-1:-2: 3]     [0, 1, 2]
+ 3 [-1:-1:-3]     []                       []                      
+ 3 [-1:-1:-3]     [0, 1, 2]
+ 3 [-1:-1:-2]     []                       []                      
+ 3 [-1:-1:-2]     [0, 1, 2]
+ 3 [-1:-1:-1]     []                       []                      
+ 3 [-1:-1:-1]     [0, 1, 2]
+ 3 [-1:-1: 1]     []                       []                      
+ 3 [-1:-1: 1]     [0, 1, 2]
+ 3 [-1:-1: 2]     []                       []                      
+ 3 [-1:-1: 2]     [0, 1, 2]
+ 3 [-1:-1: 3]     []                       []                      
+ 3 [-1:-1: 3]     [0, 1, 2]
+ 3 [-1: 0:-3]     [2]                      [2]                     
+ 3 [-1: 0:-3]     [0, 1, -1]
+ 3 [-1: 0:-2]     [2]                      [2]                     
+ 3 [-1: 0:-2]     [0, 1, -1]
+ 3 [-1: 0:-1]     [2, 1]                   [2, 1]                  
+ 3 [-1: 0:-1]     [0, -1, -1]
+ 3 [-1: 0: 1]     []                       []                      
+ 3 [-1: 0: 1]     [0, 1, 2]
+ 3 [-1: 0: 2]     []                       []                      
+ 3 [-1: 0: 2]     [0, 1, 2]
+ 3 [-1: 0: 3]     []                       []                      
+ 3 [-1: 0: 3]     [0, 1, 2]
+ 3 [-1: 1:-3]     [2]                      [2]                     
+ 3 [-1: 1:-3]     [0, 1, -1]
+ 3 [-1: 1:-2]     [2]                      [2]                     
+ 3 [-1: 1:-2]     [0, 1, -1]
+ 3 [-1: 1:-1]     [2]                      [2]                     
+ 3 [-1: 1:-1]     [0, 1, -1]
+ 3 [-1: 1: 1]     []                       []                      
+ 3 [-1: 1: 1]     [0, 1, 2]
+ 3 [-1: 1: 2]     []                       []                      
+ 3 [-1: 1: 2]     [0, 1, 2]
+ 3 [-1: 1: 3]     []                       []                      
+ 3 [-1: 1: 3]     [0, 1, 2]
+ 3 [-1: 2:-3]     []                       []                      
+ 3 [-1: 2:-3]     [0, 1, 2]
+ 3 [-1: 2:-2]     []                       []                      
+ 3 [-1: 2:-2]     [0, 1, 2]
+ 3 [-1: 2:-1]     []                       []                      
+ 3 [-1: 2:-1]     [0, 1, 2]
+ 3 [-1: 2: 1]     []                       []                      
+ 3 [-1: 2: 1]     [0, 1, 2]
+ 3 [-1: 2: 2]     []                       []                      
+ 3 [-1: 2: 2]     [0, 1, 2]
+ 3 [-1: 2: 3]     []                       []                      
+ 3 [-1: 2: 3]     [0, 1, 2]
+ 3 [-1: 3:-3]     []                       []                      
+ 3 [-1: 3:-3]     [0, 1, 2]
+ 3 [-1: 3:-2]     []                       []                      
+ 3 [-1: 3:-2]     [0, 1, 2]
+ 3 [-1: 3:-1]     []                       []                      
+ 3 [-1: 3:-1]     [0, 1, 2]
+ 3 [-1: 3: 1]     [2]                      [2]                     
+ 3 [-1: 3: 1]     [0, 1, -1]
+ 3 [-1: 3: 2]     [2]                      [2]                     
+ 3 [-1: 3: 2]     [0, 1, -1]
+ 3 [-1: 3: 3]     [2]                      [2]                     
+ 3 [-1: 3: 3]     [0, 1, -1]
+ 3 [ 0:-3:-3]     []                       []                      
+ 3 [ 0:-3:-3]     [0, 1, 2]
+ 3 [ 0:-3:-2]     []                       []                      
+ 3 [ 0:-3:-2]     [0, 1, 2]
+ 3 [ 0:-3:-1]     []                       []                      
+ 3 [ 0:-3:-1]     [0, 1, 2]
+ 3 [ 0:-3: 1]     []                       []                      
+ 3 [ 0:-3: 1]     [0, 1, 2]
+ 3 [ 0:-3: 2]     []                       []                      
+ 3 [ 0:-3: 2]     [0, 1, 2]
+ 3 [ 0:-3: 3]     []                       []                      
+ 3 [ 0:-3: 3]     [0, 1, 2]
+ 3 [ 0:-2:-3]     []                       []                      
+ 3 [ 0:-2:-3]     [0, 1, 2]
+ 3 [ 0:-2:-2]     []                       []                      
+ 3 [ 0:-2:-2]     [0, 1, 2]
+ 3 [ 0:-2:-1]     []                       []                      
+ 3 [ 0:-2:-1]     [0, 1, 2]
+ 3 [ 0:-2: 1]     [0]                      [0]                     
+ 3 [ 0:-2: 1]     [-1, 1, 2]
+ 3 [ 0:-2: 2]     [0]                      [0]                     
+ 3 [ 0:-2: 2]     [-1, 1, 2]
+ 3 [ 0:-2: 3]     [0]                      [0]                     
+ 3 [ 0:-2: 3]     [-1, 1, 2]
+ 3 [ 0:-1:-3]     []                       []                      
+ 3 [ 0:-1:-3]     [0, 1, 2]
+ 3 [ 0:-1:-2]     []                       []                      
+ 3 [ 0:-1:-2]     [0, 1, 2]
+ 3 [ 0:-1:-1]     []                       []                      
+ 3 [ 0:-1:-1]     [0, 1, 2]
+ 3 [ 0:-1: 1]     [0, 1]                   [0, 1]                  
+ 3 [ 0:-1: 1]     [-1, -1, 2]
+ 3 [ 0:-1: 2]     [0]                      [0]                     
+ 3 [ 0:-1: 2]     [-1, 1, 2]
+ 3 [ 0:-1: 3]     [0]                      [0]                     
+ 3 [ 0:-1: 3]     [-1, 1, 2]
+ 3 [ 0: 0:-3]     []                       []                      
+ 3 [ 0: 0:-3]     [0, 1, 2]
+ 3 [ 0: 0:-2]     []                       []                      
+ 3 [ 0: 0:-2]     [0, 1, 2]
+ 3 [ 0: 0:-1]     []                       []                      
+ 3 [ 0: 0:-1]     [0, 1, 2]
+ 3 [ 0: 0: 1]     []                       []                      
+ 3 [ 0: 0: 1]     [0, 1, 2]
+ 3 [ 0: 0: 2]     []                       []                      
+ 3 [ 0: 0: 2]     [0, 1, 2]
+ 3 [ 0: 0: 3]     []                       []                      
+ 3 [ 0: 0: 3]     [0, 1, 2]
+ 3 [ 0: 1:-3]     []                       []                      
+ 3 [ 0: 1:-3]     [0, 1, 2]
+ 3 [ 0: 1:-2]     []                       []                      
+ 3 [ 0: 1:-2]     [0, 1, 2]
+ 3 [ 0: 1:-1]     []                       []                      
+ 3 [ 0: 1:-1]     [0, 1, 2]
+ 3 [ 0: 1: 1]     [0]                      [0]                     
+ 3 [ 0: 1: 1]     [-1, 1, 2]
+ 3 [ 0: 1: 2]     [0]                      [0]                     
+ 3 [ 0: 1: 2]     [-1, 1, 2]
+ 3 [ 0: 1: 3]     [0]                      [0]                     
+ 3 [ 0: 1: 3]     [-1, 1, 2]
+ 3 [ 0: 2:-3]     []                       []                      
+ 3 [ 0: 2:-3]     [0, 1, 2]
+ 3 [ 0: 2:-2]     []                       []                      
+ 3 [ 0: 2:-2]     [0, 1, 2]
+ 3 [ 0: 2:-1]     []                       []                      
+ 3 [ 0: 2:-1]     [0, 1, 2]
+ 3 [ 0: 2: 1]     [0, 1]                   [0, 1]                  
+ 3 [ 0: 2: 1]     [-1, -1, 2]
+ 3 [ 0: 2: 2]     [0]                      [0]                     
+ 3 [ 0: 2: 2]     [-1, 1, 2]
+ 3 [ 0: 2: 3]     [0]                      [0]                     
+ 3 [ 0: 2: 3]     [-1, 1, 2]
+ 3 [ 0: 3:-3]     []                       []                      
+ 3 [ 0: 3:-3]     [0, 1, 2]
+ 3 [ 0: 3:-2]     []                       []                      
+ 3 [ 0: 3:-2]     [0, 1, 2]
+ 3 [ 0: 3:-1]     []                       []                      
+ 3 [ 0: 3:-1]     [0, 1, 2]
+ 3 [ 0: 3: 1]     [0, 1, 2]                [0, 1, 2]               
+ 3 [ 0: 3: 1]     [-1, -1, -1]
+ 3 [ 0: 3: 2]     [0, 2]                   [0, 2]                  
+ 3 [ 0: 3: 2]     [-1, 1, -1]
+ 3 [ 0: 3: 3]     [0]                      [0]                     
+ 3 [ 0: 3: 3]     [-1, 1, 2]
+ 3 [ 1:-3:-3]     [1]                      [1]                     
+ 3 [ 1:-3:-3]     [0, -1, 2]
+ 3 [ 1:-3:-2]     [1]                      [1]                     
+ 3 [ 1:-3:-2]     [0, -1, 2]
+ 3 [ 1:-3:-1]     [1]                      [1]                     
+ 3 [ 1:-3:-1]     [0, -1, 2]
+ 3 [ 1:-3: 1]     []                       []                      
+ 3 [ 1:-3: 1]     [0, 1, 2]
+ 3 [ 1:-3: 2]     []                       []                      
+ 3 [ 1:-3: 2]     [0, 1, 2]
+ 3 [ 1:-3: 3]     []                       []                      
+ 3 [ 1:-3: 3]     [0, 1, 2]
+ 3 [ 1:-2:-3]     []                       []                      
+ 3 [ 1:-2:-3]     [0, 1, 2]
+ 3 [ 1:-2:-2]     []                       []                      
+ 3 [ 1:-2:-2]     [0, 1, 2]
+ 3 [ 1:-2:-1]     []                       []                      
+ 3 [ 1:-2:-1]     [0, 1, 2]
+ 3 [ 1:-2: 1]     []                       []                      
+ 3 [ 1:-2: 1]     [0, 1, 2]
+ 3 [ 1:-2: 2]     []                       []                      
+ 3 [ 1:-2: 2]     [0, 1, 2]
+ 3 [ 1:-2: 3]     []                       []                      
+ 3 [ 1:-2: 3]     [0, 1, 2]
+ 3 [ 1:-1:-3]     []                       []                      
+ 3 [ 1:-1:-3]     [0, 1, 2]
+ 3 [ 1:-1:-2]     []                       []                      
+ 3 [ 1:-1:-2]     [0, 1, 2]
+ 3 [ 1:-1:-1]     []                       []                      
+ 3 [ 1:-1:-1]     [0, 1, 2]
+ 3 [ 1:-1: 1]     [1]                      [1]                     
+ 3 [ 1:-1: 1]     [0, -1, 2]
+ 3 [ 1:-1: 2]     [1]                      [1]                     
+ 3 [ 1:-1: 2]     [0, -1, 2]
+ 3 [ 1:-1: 3]     [1]                      [1]                     
+ 3 [ 1:-1: 3]     [0, -1, 2]
+ 3 [ 1: 0:-3]     [1]                      [1]                     
+ 3 [ 1: 0:-3]     [0, -1, 2]
+ 3 [ 1: 0:-2]     [1]                      [1]                     
+ 3 [ 1: 0:-2]     [0, -1, 2]
+ 3 [ 1: 0:-1]     [1]                      [1]                     
+ 3 [ 1: 0:-1]     [0, -1, 2]
+ 3 [ 1: 0: 1]     []                       []                      
+ 3 [ 1: 0: 1]     [0, 1, 2]
+ 3 [ 1: 0: 2]     []                       []                      
+ 3 [ 1: 0: 2]     [0, 1, 2]
+ 3 [ 1: 0: 3]     []                       []                      
+ 3 [ 1: 0: 3]     [0, 1, 2]
+ 3 [ 1: 1:-3]     []                       []                      
+ 3 [ 1: 1:-3]     [0, 1, 2]
+ 3 [ 1: 1:-2]     []                       []                      
+ 3 [ 1: 1:-2]     [0, 1, 2]
+ 3 [ 1: 1:-1]     []                       []                      
+ 3 [ 1: 1:-1]     [0, 1, 2]
+ 3 [ 1: 1: 1]     []                       []                      
+ 3 [ 1: 1: 1]     [0, 1, 2]
+ 3 [ 1: 1: 2]     []                       []                      
+ 3 [ 1: 1: 2]     [0, 1, 2]
+ 3 [ 1: 1: 3]     []                       []                      
+ 3 [ 1: 1: 3]     [0, 1, 2]
+ 3 [ 1: 2:-3]     []                       []                      
+ 3 [ 1: 2:-3]     [0, 1, 2]
+ 3 [ 1: 2:-2]     []                       []                      
+ 3 [ 1: 2:-2]     [0, 1, 2]
+ 3 [ 1: 2:-1]     []                       []                      
+ 3 [ 1: 2:-1]     [0, 1, 2]
+ 3 [ 1: 2: 1]     [1]                      [1]                     
+ 3 [ 1: 2: 1]     [0, -1, 2]
+ 3 [ 1: 2: 2]     [1]                      [1]                     
+ 3 [ 1: 2: 2]     [0, -1, 2]
+ 3 [ 1: 2: 3]     [1]                      [1]                     
+ 3 [ 1: 2: 3]     [0, -1, 2]
+ 3 [ 1: 3:-3]     []                       []                      
+ 3 [ 1: 3:-3]     [0, 1, 2]
+ 3 [ 1: 3:-2]     []                       []                      
+ 3 [ 1: 3:-2]     [0, 1, 2]
+ 3 [ 1: 3:-1]     []                       []                      
+ 3 [ 1: 3:-1]     [0, 1, 2]
+ 3 [ 1: 3: 1]     [1, 2]                   [1, 2]                  
+ 3 [ 1: 3: 1]     [0, -1, -1]
+ 3 [ 1: 3: 2]     [1]                      [1]                     
+ 3 [ 1: 3: 2]     [0, -1, 2]
+ 3 [ 1: 3: 3]     [1]                      [1]                     
+ 3 [ 1: 3: 3]     [0, -1, 2]
+ 3 [ 2:-3:-3]     [2]                      [2]                     
+ 3 [ 2:-3:-3]     [0, 1, -1]
+ 3 [ 2:-3:-2]     [2]                      [2]                     
+ 3 [ 2:-3:-2]     [0, 1, -1]
+ 3 [ 2:-3:-1]     [2, 1]                   [2, 1]                  
+ 3 [ 2:-3:-1]     [0, -1, -1]
+ 3 [ 2:-3: 1]     []                       []                      
+ 3 [ 2:-3: 1]     [0, 1, 2]
+ 3 [ 2:-3: 2]     []                       []                      
+ 3 [ 2:-3: 2]     [0, 1, 2]
+ 3 [ 2:-3: 3]     []                       []                      
+ 3 [ 2:-3: 3]     [0, 1, 2]
+ 3 [ 2:-2:-3]     [2]                      [2]                     
+ 3 [ 2:-2:-3]     [0, 1, -1]
+ 3 [ 2:-2:-2]     [2]                      [2]                     
+ 3 [ 2:-2:-2]     [0, 1, -1]
+ 3 [ 2:-2:-1]     [2]                      [2]                     
+ 3 [ 2:-2:-1]     [0, 1, -1]
+ 3 [ 2:-2: 1]     []                       []                      
+ 3 [ 2:-2: 1]     [0, 1, 2]
+ 3 [ 2:-2: 2]     []                       []                      
+ 3 [ 2:-2: 2]     [0, 1, 2]
+ 3 [ 2:-2: 3]     []                       []                      
+ 3 [ 2:-2: 3]     [0, 1, 2]
+ 3 [ 2:-1:-3]     []                       []                      
+ 3 [ 2:-1:-3]     [0, 1, 2]
+ 3 [ 2:-1:-2]     []                       []                      
+ 3 [ 2:-1:-2]     [0, 1, 2]
+ 3 [ 2:-1:-1]     []                       []                      
+ 3 [ 2:-1:-1]     [0, 1, 2]
+ 3 [ 2:-1: 1]     []                       []                      
+ 3 [ 2:-1: 1]     [0, 1, 2]
+ 3 [ 2:-1: 2]     []                       []                      
+ 3 [ 2:-1: 2]     [0, 1, 2]
+ 3 [ 2:-1: 3]     []                       []                      
+ 3 [ 2:-1: 3]     [0, 1, 2]
+ 3 [ 2: 0:-3]     [2]                      [2]                     
+ 3 [ 2: 0:-3]     [0, 1, -1]
+ 3 [ 2: 0:-2]     [2]                      [2]                     
+ 3 [ 2: 0:-2]     [0, 1, -1]
+ 3 [ 2: 0:-1]     [2, 1]                   [2, 1]                  
+ 3 [ 2: 0:-1]     [0, -1, -1]
+ 3 [ 2: 0: 1]     []                       []                      
+ 3 [ 2: 0: 1]     [0, 1, 2]
+ 3 [ 2: 0: 2]     []                       []                      
+ 3 [ 2: 0: 2]     [0, 1, 2]
+ 3 [ 2: 0: 3]     []                       []                      
+ 3 [ 2: 0: 3]     [0, 1, 2]
+ 3 [ 2: 1:-3]     [2]                      [2]                     
+ 3 [ 2: 1:-3]     [0, 1, -1]
+ 3 [ 2: 1:-2]     [2]                      [2]                     
+ 3 [ 2: 1:-2]     [0, 1, -1]
+ 3 [ 2: 1:-1]     [2]                      [2]                     
+ 3 [ 2: 1:-1]     [0, 1, -1]
+ 3 [ 2: 1: 1]     []                       []                      
+ 3 [ 2: 1: 1]     [0, 1, 2]
+ 3 [ 2: 1: 2]     []                       []                      
+ 3 [ 2: 1: 2]     [0, 1, 2]
+ 3 [ 2: 1: 3]     []                       []                      
+ 3 [ 2: 1: 3]     [0, 1, 2]
+ 3 [ 2: 2:-3]     []                       []                      
+ 3 [ 2: 2:-3]     [0, 1, 2]
+ 3 [ 2: 2:-2]     []                       []                      
+ 3 [ 2: 2:-2]     [0, 1, 2]
+ 3 [ 2: 2:-1]     []                       []                      
+ 3 [ 2: 2:-1]     [0, 1, 2]
+ 3 [ 2: 2: 1]     []                       []                      
+ 3 [ 2: 2: 1]     [0, 1, 2]
+ 3 [ 2: 2: 2]     []                       []                      
+ 3 [ 2: 2: 2]     [0, 1, 2]
+ 3 [ 2: 2: 3]     []                       []                      
+ 3 [ 2: 2: 3]     [0, 1, 2]
+ 3 [ 2: 3:-3]     []                       []                      
+ 3 [ 2: 3:-3]     [0, 1, 2]
+ 3 [ 2: 3:-2]     []                       []                      
+ 3 [ 2: 3:-2]     [0, 1, 2]
+ 3 [ 2: 3:-1]     []                       []                      
+ 3 [ 2: 3:-1]     [0, 1, 2]
+ 3 [ 2: 3: 1]     [2]                      [2]                     
+ 3 [ 2: 3: 1]     [0, 1, -1]
+ 3 [ 2: 3: 2]     [2]                      [2]                     
+ 3 [ 2: 3: 2]     [0, 1, -1]
+ 3 [ 2: 3: 3]     [2]                      [2]                     
+ 3 [ 2: 3: 3]     [0, 1, -1]
+ 3 [ 3:-3:-3]     [2]                      [2]                     
+ 3 [ 3:-3:-3]     [0, 1, -1]
+ 3 [ 3:-3:-2]     [2]                      [2]                     
+ 3 [ 3:-3:-2]     [0, 1, -1]
+ 3 [ 3:-3:-1]     [2, 1]                   [2, 1]                  
+ 3 [ 3:-3:-1]     [0, -1, -1]
+ 3 [ 3:-3: 1]     []                       []                      
+ 3 [ 3:-3: 1]     [0, 1, 2]
+ 3 [ 3:-3: 2]     []                       []                      
+ 3 [ 3:-3: 2]     [0, 1, 2]
+ 3 [ 3:-3: 3]     []                       []                      
+ 3 [ 3:-3: 3]     [0, 1, 2]
+ 3 [ 3:-2:-3]     [2]                      [2]                     
+ 3 [ 3:-2:-3]     [0, 1, -1]
+ 3 [ 3:-2:-2]     [2]                      [2]                     
+ 3 [ 3:-2:-2]     [0, 1, -1]
+ 3 [ 3:-2:-1]     [2]                      [2]                     
+ 3 [ 3:-2:-1]     [0, 1, -1]
+ 3 [ 3:-2: 1]     []                       []                      
+ 3 [ 3:-2: 1]     [0, 1, 2]
+ 3 [ 3:-2: 2]     []                       []                      
+ 3 [ 3:-2: 2]     [0, 1, 2]
+ 3 [ 3:-2: 3]     []                       []                      
+ 3 [ 3:-2: 3]     [0, 1, 2]
+ 3 [ 3:-1:-3]     []                       []                      
+ 3 [ 3:-1:-3]     [0, 1, 2]
+ 3 [ 3:-1:-2]     []                       []                      
+ 3 [ 3:-1:-2]     [0, 1, 2]
+ 3 [ 3:-1:-1]     []                       []                      
+ 3 [ 3:-1:-1]     [0, 1, 2]
+ 3 [ 3:-1: 1]     []                       []                      
+ 3 [ 3:-1: 1]     [0, 1, 2]
+ 3 [ 3:-1: 2]     []                       []                      
+ 3 [ 3:-1: 2]     [0, 1, 2]
+ 3 [ 3:-1: 3]     []                       []                      
+ 3 [ 3:-1: 3]     [0, 1, 2]
+ 3 [ 3: 0:-3]     [2]                      [2]                     
+ 3 [ 3: 0:-3]     [0, 1, -1]
+ 3 [ 3: 0:-2]     [2]                      [2]                     
+ 3 [ 3: 0:-2]     [0, 1, -1]
+ 3 [ 3: 0:-1]     [2, 1]                   [2, 1]                  
+ 3 [ 3: 0:-1]     [0, -1, -1]
+ 3 [ 3: 0: 1]     []                       []                      
+ 3 [ 3: 0: 1]     [0, 1, 2]
+ 3 [ 3: 0: 2]     []                       []                      
+ 3 [ 3: 0: 2]     [0, 1, 2]
+ 3 [ 3: 0: 3]     []                       []                      
+ 3 [ 3: 0: 3]     [0, 1, 2]
+ 3 [ 3: 1:-3]     [2]                      [2]                     
+ 3 [ 3: 1:-3]     [0, 1, -1]
+ 3 [ 3: 1:-2]     [2]                      [2]                     
+ 3 [ 3: 1:-2]     [0, 1, -1]
+ 3 [ 3: 1:-1]     [2]                      [2]                     
+ 3 [ 3: 1:-1]     [0, 1, -1]
+ 3 [ 3: 1: 1]     []                       []                      
+ 3 [ 3: 1: 1]     [0, 1, 2]
+ 3 [ 3: 1: 2]     []                       []                      
+ 3 [ 3: 1: 2]     [0, 1, 2]
+ 3 [ 3: 1: 3]     []                       []                      
+ 3 [ 3: 1: 3]     [0, 1, 2]
+ 3 [ 3: 2:-3]     []                       []                      
+ 3 [ 3: 2:-3]     [0, 1, 2]
+ 3 [ 3: 2:-2]     []                       []                      
+ 3 [ 3: 2:-2]     [0, 1, 2]
+ 3 [ 3: 2:-1]     []                       []                      
+ 3 [ 3: 2:-1]     [0, 1, 2]
+ 3 [ 3: 2: 1]     []                       []                      
+ 3 [ 3: 2: 1]     [0, 1, 2]
+ 3 [ 3: 2: 2]     []                       []                      
+ 3 [ 3: 2: 2]     [0, 1, 2]
+ 3 [ 3: 2: 3]     []                       []                      
+ 3 [ 3: 2: 3]     [0, 1, 2]
+ 3 [ 3: 3:-3]     []                       []                      
+ 3 [ 3: 3:-3]     [0, 1, 2]
+ 3 [ 3: 3:-2]     []                       []                      
+ 3 [ 3: 3:-2]     [0, 1, 2]
+ 3 [ 3: 3:-1]     []                       []                      
+ 3 [ 3: 3:-1]     [0, 1, 2]
+ 3 [ 3: 3: 1]     []                       []                      
+ 3 [ 3: 3: 1]     [0, 1, 2]
+ 3 [ 3: 3: 2]     []                       []                      
+ 3 [ 3: 3: 2]     [0, 1, 2]
+ 3 [ 3: 3: 3]     []                       []                      
+ 3 [ 3: 3: 3]     [0, 1, 2]

tests/slicing2.py

@@ -0,0 +1,8 @@
+try:
+    import ulab as np
+except:
+    import numpy as np
+
+a = np.array(range(9), dtype=np.float)
+print("a:\t", list(a))
+print("a < 5:\t", list(a[a < 5]))

tests/slicing2.py.exp

@@ -0,0 +1,2 @@
+a:	 [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+a < 5:	 [0.0, 1.0, 2.0, 3.0, 4.0]