Testing

README.md

@@ -8,9 +8,9 @@ Documentation can be found under https://micropython-ulab.readthedocs.io/en/late
 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.
 
-Firmware for pyboard.v.1.1 is updated once in a while, and can be downloaded
-from https://github.com/v923z/micropython-ulab/releases, otherwise, it can 
-be compiled from the source by following the steps
+Firmware for pyboard.v.1.1, as well as for PYBD_SF6 is updated once 
+in a while, and can be downloaded from https://github.com/v923z/micropython-ulab/releases, 
+otherwise, it can be compiled from the source by following the steps
 https://micropython-usermod.readthedocs.io/en/latest/usermods_05.html#compiling-our-module.
 
 

code/fft.c

@@ -26,13 +26,13 @@ enum FFT_TYPE {
     FFT_SPECTRUM,
 };
 
-void fft_kernel(float *real, float *imag, int n, int isign) {
+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. 
     int j, m, mmax, istep;
-    float tempr, tempi;
-    float wtemp, wr, wpr, wpi, wi, theta;
+    mp_float_t tempr, tempi;
+    mp_float_t wtemp, wr, wpr, wpi, wi, theta;
 
     j = 0;
     for(int i = 0; i < n; i++) {
@@ -52,9 +52,9 @@ void fft_kernel(float *real, float *imag, int n, int isign) {
     while (n > mmax) {
         istep = mmax << 1;
         theta = -1.0*isign*6.28318530717959/istep;
-        wtemp = sinf(0.5 * theta);
+        wtemp = MICROPY_FLOAT_C_FUN(sin)(0.5 * theta);
         wpr = -2.0 * wtemp * wtemp;
-        wpi = sinf(theta);
+        wpi = MICROPY_FLOAT_C_FUN(sin)(theta);
         wr = 1.0;
         wi = 0.0;
         for(m = 0; m < mmax; m++) {
@@ -92,27 +92,27 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
     }
 
     ndarray_obj_t *out_re = create_new_ndarray(1, len, NDARRAY_FLOAT);
-    float *data_re = (float *)out_re->array->items;
+    mp_float_t *data_re = (mp_float_t *)out_re->array->items;
 
     if(re->array->typecode == NDARRAY_FLOAT) { 
         // By treating this case separately, we can save a bit of time.
         // I don't know if it is worthwhile, though...
-        memcpy((float *)out_re->array->items, (float *)re->array->items, re->bytes);
+        memcpy((mp_float_t *)out_re->array->items, (mp_float_t *)re->array->items, re->bytes);
     } else {
         for(size_t i=0; i < len; i++) {
             data_re[i] = ndarray_get_float_value(re->array->items, re->array->typecode, i);
         }
     }
     ndarray_obj_t *out_im = create_new_ndarray(1, len, NDARRAY_FLOAT);
-    float *data_im = (float *)out_im->array->items;
+    mp_float_t *data_im = (mp_float_t *)out_im->array->items;
 
     if(n_args == 2) {
         ndarray_obj_t *im = MP_OBJ_TO_PTR(arg_im);
         if (re->array->len != im->array->len) {
             mp_raise_ValueError("real and imaginary parts must be of equal length");
         }
         if(im->array->typecode == NDARRAY_FLOAT) {
-            memcpy((float *)out_im->array->items, (float *)im->array->items, im->bytes);
+            memcpy((mp_float_t *)out_im->array->items, (mp_float_t *)im->array->items, im->bytes);
         } else {
             for(size_t i=0; i < len; i++) {
                 data_im[i] = ndarray_get_float_value(im->array->items, im->array->typecode, i);
@@ -123,7 +123,7 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
         fft_kernel(data_re, data_im, len, 1);
         if(type == FFT_SPECTRUM) {
             for(size_t i=0; i < len; i++) {
-                data_re[i] = sqrtf(data_re[i]*data_re[i] + data_im[i]*data_im[i]);
+                data_re[i] = MICROPY_FLOAT_C_FUN(sqrt)(data_re[i]*data_re[i] + data_im[i]*data_im[i]);
             }
         }
     } else { // inverse transform

code/linalg.c

@@ -108,24 +108,24 @@ mp_obj_t linalg_size(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args)
     }
 }
 
-bool linalg_invert_matrix(float *data, size_t N) {
+bool linalg_invert_matrix(mp_float_t *data, size_t N) {
     // returns true, of the inversion was successful, 
     // false, if the matrix is singular
 
     // initially, this is the unit matrix: the contents of this matrix is what 
     // will be returned after all the transformations
-    float *unit = m_new(float, N*N);
+    mp_float_t *unit = m_new(mp_float_t, N*N);
 
-    float elem = 1.0;
+    mp_float_t elem = 1.0;
     // initialise the unit matrix
-    memset(unit, 0, sizeof(float)*N*N);
+    memset(unit, 0, sizeof(mp_float_t)*N*N);
     for(size_t m=0; m < N; m++) {
-        memcpy(&unit[m*(N+1)], &elem, sizeof(float));
+        memcpy(&unit[m*(N+1)], &elem, sizeof(mp_float_t));
     }
     for(size_t m=0; m < N; m++){
         // this could be faster with ((c < epsilon) && (c > -epsilon))
         if(abs(data[m*(N+1)]) < epsilon) {
-            m_del(float, unit, N*N);
+            m_del(mp_float_t, unit, N*N);
             return false;
         }
         for(size_t n=0; n < N; n++){
@@ -145,8 +145,8 @@ bool linalg_invert_matrix(float *data, size_t N) {
             unit[N*m+n] /= elem;
         }
     }
-    memcpy(data, unit, sizeof(float)*N*N);
-    m_del(float, unit, N*N);
+    memcpy(data, unit, sizeof(mp_float_t)*N*N);
+    m_del(mp_float_t, unit, N*N);
     return true;
 }
 
@@ -163,21 +163,21 @@ mp_obj_t linalg_inv(mp_obj_t o_in) {
         mp_raise_ValueError("only square matrices can be inverted");
     }
     ndarray_obj_t *inverted = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT);
-    float *data = (float *)inverted->array->items;
+    mp_float_t *data = (mp_float_t *)inverted->array->items;
     mp_obj_t elem;
     for(size_t m=0; m < o->m; m++) { // rows first
         for(size_t n=0; n < o->n; n++) { // columns next
             // this could, perhaps, be done in single line... 
             // On the other hand, we probably spend little time here
             elem = mp_binary_get_val_array(o->array->typecode, o->array->items, m*o->n+n);
-            data[m*o->n+n] = (float)mp_obj_get_float(elem);
+            data[m*o->n+n] = (mp_float_t)mp_obj_get_float(elem);
         }
     }
 
     if(!linalg_invert_matrix(data, o->m)) {
         // TODO: I am not sure this is needed here. Otherwise, 
         // how should we free up the unused RAM of inverted?
-        m_del(float, inverted->array->items, o->n*o->n);
+        m_del(mp_float_t, inverted->array->items, o->n*o->n);
         mp_raise_ValueError("input matrix is singular");
     }
     return MP_OBJ_FROM_PTR(inverted);
@@ -192,8 +192,8 @@ 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);
-    float *outdata = (float *)out->array->items;
-    float sum, v1, v2;
+    mp_float_t *outdata = (mp_float_t *)out->array->items;
+    mp_float_t sum, v1, v2;
     for(size_t i=0; i < m1->n; i++) {
         for(size_t j=0; j < m2->m; j++) {
             sum = 0.0;
@@ -301,14 +301,14 @@ mp_obj_t linalg_det(mp_obj_t oin) {
         mp_raise_ValueError("input must be square matrix");
     }
 
-    float *tmp = m_new(float, in->n*in->n);
+    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);
     }
-    float c;
+    mp_float_t c;
     for(size_t m=0; m < in->m-1; m++){
         if(abs(tmp[m*(in->n+1)]) < epsilon) {
-            m_del(float, tmp, in->n*in->n);
+            m_del(mp_float_t, tmp, in->n*in->n);
             return mp_obj_new_float(0.0);
         }
         for(size_t n=0; n < in->n; n++){
@@ -320,12 +320,12 @@ mp_obj_t linalg_det(mp_obj_t oin) {
             }
         }
     }
-    float det = 1.0;
+    mp_float_t det = 1.0;
 
     for(size_t m=0; m < in->m; m++){ 
         det *= tmp[m*(in->n+1)];
     }
-    m_del(float, tmp, in->n*in->n);
+    m_del(mp_float_t, tmp, in->n*in->n);
     return mp_obj_new_float(det);
 }
 
@@ -337,14 +337,15 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
     if(in->m != in->n) {
         mp_raise_ValueError("input must be square matrix");
     }
-    float *array = m_new(float, in->array->len);
+    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);
     }
     // make sure the matrix is symmetric
     for(size_t m=0; m < in->m; m++) {
         for(size_t n=m+1; n < in->n; n++) {
             // compare entry (m, n) to (n, m)
+            // TODO: this must probably be scaled!
             if(epsilon < abs(array[m*in->n + n] - array[n*in->n + m])) {
                 mp_raise_ValueError("input matrix is asymmetric");
             }
@@ -354,12 +355,12 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
     // if we got this far, then the matrix will be symmetric
 
     ndarray_obj_t *eigenvectors = create_new_ndarray(in->m, in->n, NDARRAY_FLOAT);
-    float *eigvectors = (float *)eigenvectors->array->items;
+    mp_float_t *eigvectors = (mp_float_t *)eigenvectors->array->items;
     // start out with the unit matrix
     for(size_t m=0; m < in->m; m++) {
         eigvectors[m*(in->n+1)] = 1.0;
     }
-    float largest, w, t, c, s, tau, aMk, aNk, vm, vn;
+    mp_float_t largest, w, t, c, s, tau, aMk, aNk, vm, vn;
     size_t M, N;
     size_t iterations = JACOBI_MAX*in->n*in->n;
     do {
@@ -387,12 +388,12 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
         // The following if/else chooses the smaller absolute value for the tangent 
         // of the rotation angle. Going with the smaller should be numerically stabler.
         if(w > 0) {
-            t = sqrtf(w*w + 1.0) - w;
+            t = MICROPY_FLOAT_C_FUN(sqrt)(w*w + 1.0) - w;
         } else {
-            t = -1.0*(sqrtf(w*w + 1.0) + w);
+            t = -1.0*(MICROPY_FLOAT_C_FUN(sqrt)(w*w + 1.0) + w);
         }
-        s = t / sqrtf(t*t + 1.0); // the sine of the rotation angle
-        c = 1.0 / sqrtf(t*t + 1.0); // the cosine of the rotation angle
+        s = t / MICROPY_FLOAT_C_FUN(sqrt)(t*t + 1.0); // the sine of the rotation angle
+        c = 1.0 / MICROPY_FLOAT_C_FUN(sqrt)(t*t + 1.0); // the cosine of the rotation angle
         tau = (1.0-c)/s; // this is equal to the tangent of the half of the rotation angle
 
         // at this point, we have the rotation angles, so we can transform the matrix
@@ -438,15 +439,15 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
 
     if(iterations == 0) { 
         // the computation did not converge; numpy raises LinAlgError
-        m_del(float, array, in->array->len);
+        m_del(mp_float_t, array, in->array->len);
         mp_raise_ValueError("iterations did not converge");
     }
     ndarray_obj_t *eigenvalues = create_new_ndarray(1, in->n, NDARRAY_FLOAT);
-    float *eigvalues = (float *)eigenvalues->array->items;
+    mp_float_t *eigvalues = (mp_float_t *)eigenvalues->array->items;
     for(size_t i=0; i < in->n; i++) {
         eigvalues[i] = array[i*(in->n+1)];
     }
-    m_del(float, array, in->array->len);
+    m_del(mp_float_t, array, in->array->len);
 
     mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(mp_obj_new_tuple(2, NULL));
     tuple->items[0] = MP_OBJ_FROM_PTR(eigenvalues);

code/linalg.h

@@ -14,13 +14,19 @@
 #include "ndarray.h"
 
 #define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }
-#define epsilon        1e-6
+
+#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
+#define epsilon        1.2e-7
+#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
+#define epsilon        2.3e-16
+#endif
+
 #define JACOBI_MAX     20
 
 mp_obj_t linalg_transpose(mp_obj_t );
 mp_obj_t linalg_reshape(mp_obj_t , mp_obj_t );
 mp_obj_t linalg_size(size_t , const mp_obj_t *, mp_map_t *);
-bool linalg_invert_matrix(float *, size_t );
+bool linalg_invert_matrix(mp_float_t *, size_t );
 mp_obj_t linalg_inv(mp_obj_t );
 mp_obj_t linalg_dot(mp_obj_t , mp_obj_t );
 mp_obj_t linalg_zeros(size_t , const mp_obj_t *, mp_map_t *);