From 2c4811cea42daae6b40a49e042a33b746075e8c2 Mon Sep 17 00:00:00 2001
From: feng-arch <feng-arch@outlook.com>
Date: Mon, 13 Oct 2025 15:42:15 +0800
Subject: [PATCH] feat: add support for eigenmath

---
 ports/rp2/boards/RPI_PICO/mpconfigboard.h |    2 +-
 user_modules/eigenmath/eheap.c            |  341 ++++
 user_modules/eigenmath/eheap.h            |   33 +
 user_modules/eigenmath/eigenmath.c        | 1932 +++++++++------------
 user_modules/eigenmath/eigenmath.h        |   35 +
 user_modules/eigenmath/eigenmath_mpy.c    |  281 +++
 user_modules/eigenmath/micropython.cmake  |    6 +-
 user_modules/eigenmath/micropython.mk     |    8 +-
 8 files changed, 1490 insertions(+), 1148 deletions(-)
 create mode 100644 user_modules/eigenmath/eheap.c
 create mode 100644 user_modules/eigenmath/eheap.h
 create mode 100644 user_modules/eigenmath/eigenmath.h
 create mode 100644 user_modules/eigenmath/eigenmath_mpy.c

diff --git a/ports/rp2/boards/RPI_PICO/mpconfigboard.h b/ports/rp2/boards/RPI_PICO/mpconfigboard.h
index 7f34cffce..43d167b88 100644
--- a/ports/rp2/boards/RPI_PICO/mpconfigboard.h
+++ b/ports/rp2/boards/RPI_PICO/mpconfigboard.h
@@ -2,4 +2,4 @@
 #define MICROPY_HW_BOARD_NAME                   "Raspberry Pi Pico"
 // Modified from MPY origin to reduce flash storage to accommodate larger program flash requirement
 // of lvgl and its bindings. Developers should review this setting when adding additional features 
-#define MICROPY_HW_FLASH_STORAGE_BYTES          (1024 * 1024)
+#define MICROPY_HW_FLASH_STORAGE_BYTES          (512 * 1024)
diff --git a/user_modules/eigenmath/eheap.c b/user_modules/eigenmath/eheap.c
new file mode 100644
index 000000000..d61564b13
--- /dev/null
+++ b/user_modules/eigenmath/eheap.c
@@ -0,0 +1,341 @@
+/*======================================================================== 
+ *  eheap.c – Improved standalone heap manager 
+ *            32-bit, aligned, safe operations
+ *======================================================================*/
+
+ #include "eheap.h"
+ #include <stdint.h>
+ #include <string.h>
+ #include <stdbool.h>
+ #include "py/runtime.h"
+ 
+ #ifndef EHEAP_ALIGN
+ #define EHEAP_ALIGN 4               /* default byte alignment (power of two) */
+ #endif
+ 
+ #define ALIGN_MASK    (EHEAP_ALIGN - 1)
+ #define ALIGN_UP(x)   (((x) + ALIGN_MASK) & ~ALIGN_MASK)
+ 
+ /*--------------------------------------------------  block header layout */
+ typedef struct block_link {
+     struct block_link *next;
+     size_t             size;     /* MSB=1 => allocated, lower bits => block size */
+ } block_t;
+ 
+ #define USED_MASK     ((size_t)1 << (sizeof(size_t)*8 - 1))
+ #define IS_USED(b)    (((b)->size) &  USED_MASK)
+ #define MARK_USED(b)  ((b)->size |= USED_MASK)
+ #define MARK_FREE(b)  ((b)->size &= ~USED_MASK)
+ #define BLOCK_SIZE(b) ((b)->size & ~USED_MASK)
+ 
+ #define HDR_SIZE     ALIGN_UP(sizeof(block_t))
+ #define MIN_SPLIT    (HDR_SIZE * 2)
+ 
+ /*--------------------------------------------------  heap globals */
+ static uint8_t *heap_base  = NULL;
+ static uint8_t *heap_end   = NULL;
+ static size_t   heap_total = 0;
+ 
+ static block_t  start_node;         /* dummy head */
+ static block_t  end_marker;         /* tail sentinel storage */
+ static block_t *end_node = &end_marker;
+ 
+ static size_t free_bytes = 0;
+ static size_t min_free   = 0;
+ static bool   initialized = false;
+ 
+ /*-------------------------------------------------------------------*/
+ static bool is_valid_block(block_t *blk) {
+     uint8_t *ptr = (uint8_t*)blk;
+     if (ptr < heap_base || ptr >= heap_end) return false;
+     return (((uintptr_t)ptr - (uintptr_t)heap_base) & ALIGN_MASK) == 0;
+ }
+ 
+/* Insert and coalesce a free block (address-ordered), with overflow guards */  
+static void insert_free(block_t *blk) {  
+    if (!is_valid_block(blk) || IS_USED(blk)) {  
+        mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("insert_free: invalid or used block"));  
+        return;  
+    }  
+    size_t blk_sz = BLOCK_SIZE(blk);  
+    /* guard pointer addition overflow */  
+    if (blk_sz > (size_t)(heap_end - (uint8_t*)blk)) {  
+        mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("insert_free: block size overflow"));  
+        return;  
+    }  
+    uint8_t *blk_end = (uint8_t*)blk + blk_sz;  
+  
+    block_t *prev = &start_node;  
+    /* find insertion point */  
+    while (prev->next < blk && prev->next != end_node) {  
+        prev = prev->next;  
+        if (!is_valid_block(prev)) {  
+            mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("insert_free: corrupted free list"));  
+            return;  
+        }  
+    }  
+  
+    /* forward merge */  
+    if (prev->next != end_node) {  
+        block_t *fwd = prev->next;  
+        if (!IS_USED(fwd) &&  
+            (uint8_t*)fwd == blk_end) {  
+            /* fuse sizes */  
+            size_t total = blk_sz + BLOCK_SIZE(fwd);  
+            if (total < blk_sz) { /* overflow? */  
+                mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("insert_free: combine overflow"));  
+                return;  
+            }  
+            blk->size = total;  
+            blk->next = fwd->next;  
+            blk_sz = total;  /* update for potential backward merge */  
+        } else {  
+            blk->next = fwd;  
+        }  
+    } else {  
+        blk->next = end_node;  
+    }  
+  
+    /* backward merge */  
+    if (prev != &start_node && !IS_USED(prev)) {  
+        uint8_t *prev_end = (uint8_t*)prev + BLOCK_SIZE(prev);  
+        if (prev_end == (uint8_t*)blk) {  
+            size_t total = BLOCK_SIZE(prev) + blk_sz;  
+            if (total < BLOCK_SIZE(prev)) { /* overflow? */  
+                mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("insert_free: combine overflow"));  
+                return;  
+            }  
+            prev->size = total;  
+            prev->next = blk->next;  
+            return;  
+        }  
+    }  
+    prev->next = blk;  
+}  
+ 
+ static void heap_init_once(void) {
+     if (initialized) return;
+ 
+     /* single free block covers [heap_base .. heap_end) excluding end_node */
+     block_t *first = (block_t*)heap_base;
+     first->size = (heap_total - HDR_SIZE);
+     MARK_FREE(first);
+     first->next = end_node;
+ 
+     start_node.next = first;
+     start_node.size = 0;
+ 
+     /* initialize end marker */
+     end_node->next = NULL;
+     end_node->size = 0;
+     MARK_USED(end_node);
+ 
+     free_bytes = BLOCK_SIZE(first);
+     min_free   = free_bytes;
+     initialized = true;
+ }
+ 
+ void eheap_init(void *buf, size_t bytes) {
+     if (!buf || bytes <= HDR_SIZE*2 + ALIGN_MASK) {
+         mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("eheap_init: invalid region"));
+         return;
+     }
+ 
+     /* align base upward */
+     uintptr_t start = ALIGN_UP((uintptr_t)buf);
+     size_t   loss  = start - (uintptr_t)buf;
+     bytes = (bytes > loss) ? bytes - loss : 0;
+     bytes = (bytes / EHEAP_ALIGN) * EHEAP_ALIGN;
+     if (bytes <= HDR_SIZE*2) {
+         mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("eheap_init: too small after align"));
+         return;
+     }
+ 
+     heap_base  = (uint8_t*)start;
+     heap_total = bytes;
+ 
+     /* reserve tail for end_node */
+     size_t res = ALIGN_UP(sizeof(block_t));
+     heap_total -= res;
+     end_node = (block_t*)(heap_base + heap_total);
+ 
+     heap_end = heap_base + heap_total;
+ 
+     initialized = false;
+     heap_init_once();
+ }
+ 
+ void* e_malloc(size_t size) {
+     if (size == 0 || !initialized) return NULL;
+ 
+     /* check overflow */
+     if (size > SIZE_MAX - HDR_SIZE) return NULL;
+     size_t needed = ALIGN_UP(size + HDR_SIZE);
+ 
+     block_t *prev = &start_node;
+     block_t *cur  = start_node.next;
+ 
+     while (cur != end_node) {
+         if (!is_valid_block(cur)) {
+             mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("e_malloc: corrupted heap"));
+             return NULL;
+         }
+         if (!IS_USED(cur) && BLOCK_SIZE(cur) >= needed) {
+             size_t remain = BLOCK_SIZE(cur) - needed;
+             if (remain >= MIN_SPLIT) {
+                 /* split */
+                 block_t *split = (block_t*)((uint8_t*)cur + needed);
+                 split->size = remain;
+                 MARK_FREE(split);
+                 split->next = cur->next;
+ 
+                 cur->size = needed;
+                 prev->next = split;
+             } else {
+                 /* use entire */
+                 prev->next = cur->next;
+                 needed = BLOCK_SIZE(cur);
+             }
+             MARK_USED(cur);
+ 
+             free_bytes -= needed;
+             if (free_bytes < min_free) min_free = free_bytes;
+             return (uint8_t*)cur + HDR_SIZE;
+         }
+         prev = cur;
+         cur  = cur->next;
+     }
+     return NULL;
+ }
+ 
+ void e_free(void *ptr) {
+     if (!ptr || !initialized) return;
+ 
+     uint8_t *p = (uint8_t*)ptr;
+     if (p < heap_base + HDR_SIZE || p >= heap_end) {
+         mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("e_free: invalid ptr"));
+         return;
+     }
+     if (((uintptr_t)p - HDR_SIZE) & ALIGN_MASK) {
+         mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("e_free: unaligned ptr"));
+         return;
+     }
+ 
+     block_t *blk = (block_t*)(p - HDR_SIZE);
+     if (!IS_USED(blk)) {
+         //mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("e_free: double free"));
+         return;
+     }
+ 
+     size_t sz = BLOCK_SIZE(blk);
+     if (sz == 0 || (uint8_t*)blk + sz > heap_end) {
+         mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("e_free: bad size"));
+         return;
+     }
+ 
+     MARK_FREE(blk);
+     free_bytes += sz;
+     insert_free(blk);
+ }
+ 
+ void* e_realloc(void *ptr, size_t new_size) {
+     if (!ptr) return e_malloc(new_size);
+     if (new_size == 0) { e_free(ptr); return NULL; }
+     if (!initialized) return NULL;
+ 
+     uint8_t *p = (uint8_t*)ptr;
+     if (p < heap_base + HDR_SIZE || p >= heap_end) {
+         mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("e_realloc: invalid ptr"));
+         return NULL;
+     }
+ 
+     block_t *blk = (block_t*)(p - HDR_SIZE);
+     if (!IS_USED(blk)) {
+         mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("e_realloc: block not used"));
+         return NULL;
+     }
+ 
+     size_t curr = BLOCK_SIZE(blk) - HDR_SIZE;
+     if (new_size <= curr) return ptr;
+ 
+     /* try expand into next free block */
+     uint8_t *next_addr = (uint8_t*)blk + BLOCK_SIZE(blk);
+     if (next_addr + HDR_SIZE <= heap_end) {
+         block_t *next = (block_t*)next_addr;
+         if (is_valid_block(next) && !IS_USED(next)) {
+             size_t combined = BLOCK_SIZE(blk) + BLOCK_SIZE(next);
+             size_t need = ALIGN_UP(new_size + HDR_SIZE);
+             if (combined >= need) {
+                 /* remove next from free list */
+                 block_t *prev = &start_node;
+                 while (prev->next != next && prev->next != end_node) {
+                     prev = prev->next;
+                 }
+                 if (prev->next == next) {
+                     prev->next = next->next;
+                     /* compute new free_bytes: remove next size */
+                     free_bytes -= BLOCK_SIZE(next);
+ 
+                     /* update blk size */
+                     blk->size = (blk->size & USED_MASK) | need;
+                     size_t leftover = combined - need;
+                     if (leftover >= MIN_SPLIT) {
+                         block_t *split = (block_t*)((uint8_t*)blk + need);
+                         split->size = leftover;
+                         MARK_FREE(split);
+                         insert_free(split);
+                     } else {
+                         /* absorb all */
+                         blk->size = (blk->size & USED_MASK) | combined;
+                     }
+                     if (free_bytes < min_free) min_free = free_bytes;
+                     return ptr;
+                 }
+             }
+         }
+     }
+ 
+     /* fallback: alloc-copy-free */
+     void *nptr = e_malloc(new_size);
+     if (nptr) {
+         memcpy(nptr, ptr, curr);
+         e_free(ptr);
+     }
+     return nptr;
+ }
+ 
+ size_t e_heap_free(void)     { return free_bytes; }
+ size_t e_heap_min_free(void) { return min_free;   }
+ 
+
+ int e_heap_fragmentation(void) {
+     if (!initialized || free_bytes == 0) return 0;
+     size_t largest = 0;
+     for (block_t *b = start_node.next; b != end_node; b = b->next) {
+         if (!IS_USED(b) && BLOCK_SIZE(b) > largest) {
+             largest = BLOCK_SIZE(b);
+         }
+     }
+     if (largest == 0) return 100;
+     return (int)(100 - (largest * 100) / free_bytes);
+ }
+ 
+ bool e_heap_validate(void) {
+     if (!initialized) return false;
+     size_t counted = 0;
+     for (block_t *b = start_node.next; b != end_node; b = b->next) {
+         if (!is_valid_block(b)) return false;
+         if ((uint8_t*)b + BLOCK_SIZE(b) > heap_end) return false;
+         if (!IS_USED(b)) {
+             counted += BLOCK_SIZE(b);
+             /* ensure no adjacent free blocks */
+             block_t *n = b->next;
+             if (n != end_node && !IS_USED(n) &&
+                 (uint8_t*)b + BLOCK_SIZE(b) == (uint8_t*)n) {
+                 return false;
+             }
+         }
+     }
+     return (counted == free_bytes);
+ }
+ 
\ No newline at end of file
diff --git a/user_modules/eigenmath/eheap.h b/user_modules/eigenmath/eheap.h
new file mode 100644
index 000000000..09545cd38
--- /dev/null
+++ b/user_modules/eigenmath/eheap.h
@@ -0,0 +1,33 @@
+/*========================================================================
+ *  eheap.h – Minimal standalone heap manager (Freertos heap_4 style)
+ *  ---------------------------------------------------------------------
+ *  API:
+ *      void   eheap_init(void *buffer, size_t size);
+ *      void*  e_malloc(size_t bytes);
+ *      void   e_free(void *ptr);
+ *      void*  e_realloc(void *ptr, size_t new_size);
+ *      size_t e_heap_free(void);
+ *      size_t e_heap_min_free(void);
+ *======================================================================*/
+
+ #ifndef EHEAP_H
+ #define EHEAP_H
+ 
+ #include <stddef.h>
+ 
+ #ifdef __cplusplus
+ extern "C" {
+ #endif
+ 
+ void   eheap_init(void *buffer, size_t size);
+ void*  e_malloc(size_t size);
+ void   e_free(void *ptr);
+ void*  e_realloc(void *ptr, size_t new_size);
+ size_t e_heap_free(void);
+ size_t e_heap_min_free(void);
+ int e_heap_fragmentation(void);
+ #ifdef __cplusplus
+ }
+ #endif
+ 
+ #endif /* EHEAP_H */
\ No newline at end of file
diff --git a/user_modules/eigenmath/eigenmath.c b/user_modules/eigenmath/eigenmath.c
index 9d8ca8039..80c6c593b 100644
--- a/user_modules/eigenmath/eigenmath.c
+++ b/user_modules/eigenmath/eigenmath.c
@@ -1,12 +1,3 @@
-// Include MicroPython API.
-#include "py/runtime.h"
-
-// Used to get the time in the Timer class example.
-#include "py/mphal.h"
-#include "py/misc.h"
-#include "shared/readline/readline.h"
-
-
 /*
 BSD 2-Clause License
 
@@ -34,6 +25,7 @@ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
+#include <stdbool.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdint.h>
@@ -44,35 +36,25 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <setjmp.h>
 #include <math.h>
 #include <errno.h>
-#include <time.h>
+#include "eigenmath.h"
+#include "eheap.h"
+#include "py/obj.h"
+//#include "py/mpconfig.h"
+#include "py/misc.h"
+#include "py/runtime.h"
+#include "py/objstr.h"
+#include "eigenmath.h"
+//#define STACKSIZE 100000 // evaluation stack
+//#define BLOCKSIZE 10000
+//#define MAXBLOCKS 2000
+//#define BUCKETSIZE 100
+//#define STRBUFLEN 1000
+//#define MAXDIM 24
 
-// #define STACKSIZE 100000 // evaluation stack
-// #define BLOCKSIZE 10000
-// #define MAXBLOCKS 2000
-// #define BUCKETSIZE 100
-// #define STRBUFLEN 1000
-// #define MAXDIM 24
 
-// #define STACKSIZE 4096    // 栈大小从100000减到4096
-// #define BLOCKSIZE 512     // 块大小从10000减到512
-// #define MAXBLOCKS 64      // 最大块数从2000减到64
-// #define BUCKETSIZE 32     // 哈希桶大小相应减小
-// #define STRBUFLEN 256     // 字符串缓冲区减小
-// #define MAXDIM 16         // 维度限制适当降低
+uint32_t STACKSIZE ;
+uint32_t MAXATOMS ; // 10,240 atoms
 
-// #define STACKSIZE 128       // 栈大小：128个元素（每个元素按4字节算约512字节）
-// #define BLOCKSIZE 16        // 块大小：每个块包含16个atom
-// #define MAXBLOCKS 64        // 最大块数：64块（16*64=1024个atom）
-// #define BUCKETSIZE 4        // 哈希桶大小：仅保留4个桶
-// #define STRBUFLEN 32        // 字符串缓冲区：每个字符串最大32字节
-// #define MAXDIM 4            // 维度限制：降低到4维
-
-#define STACKSIZE 512     // 进一步减小栈大小
-#define BLOCKSIZE 128     // 减小块大小  
-#define MAXBLOCKS 64      // 减少最大块数
-#define BUCKETSIZE 20     // 减小哈希桶
-#define STRBUFLEN 256      // 减小字符串缓冲区
-#define MAXDIM 16          // 降低维度限制
 
 // MAXBLOCKS * BLOCKSIZE = 20,000,000 atoms
 
@@ -101,7 +83,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 //              |SYM    |    |SYM    |    |SYM    |
 //              |"mul"  |    |"a"    |    |"b"    |
 //              |_______|    |_______|    |_______|
-
+/**/
 struct atom {
 	union {
 		struct {
@@ -160,7 +142,6 @@ struct tensor {
 #define ARG		(0 * BUCKETSIZE + 9)
 
 #define BINDING		(1 * BUCKETSIZE + 0)
-#define BREAK		(1 * BUCKETSIZE + 1)
 
 #define C_UPPER		(2 * BUCKETSIZE + 0)
 #define C_LOWER		(2 * BUCKETSIZE + 1)
@@ -223,9 +204,7 @@ struct tensor {
 #define KRONECKER	(10 * BUCKETSIZE + 0)
 
 #define LAST		(11 * BUCKETSIZE + 0)
-#define LGAMMA		(11 * BUCKETSIZE + 1)
-#define LOG		(11 * BUCKETSIZE + 2)
-#define LOOP		(11 * BUCKETSIZE + 3)
+#define LOG		(11 * BUCKETSIZE + 1)
 
 #define MAG		(12 * BUCKETSIZE + 0)
 #define MINOR		(12 * BUCKETSIZE + 1)
@@ -256,14 +235,13 @@ struct tensor {
 
 #define R_UPPER		(17 * BUCKETSIZE + 0)
 #define R_LOWER		(17 * BUCKETSIZE + 1)
-#define RAND		(17 * BUCKETSIZE + 2)
-#define RANK		(17 * BUCKETSIZE + 3)
-#define RATIONALIZE	(17 * BUCKETSIZE + 4)
-#define REAL		(17 * BUCKETSIZE + 5)
-#define RECTF		(17 * BUCKETSIZE + 6)
-#define ROOTS		(17 * BUCKETSIZE + 7)
-#define ROTATE		(17 * BUCKETSIZE + 8)
-#define RUN		(17 * BUCKETSIZE + 9)
+#define RANK		(17 * BUCKETSIZE + 2)
+#define RATIONALIZE	(17 * BUCKETSIZE + 3)
+#define REAL		(17 * BUCKETSIZE + 4)
+#define RECTF		(17 * BUCKETSIZE + 5)
+#define ROOTS		(17 * BUCKETSIZE + 6)
+#define ROTATE		(17 * BUCKETSIZE + 7)
+#define RUN		(17 * BUCKETSIZE + 8)
 
 #define S_UPPER		(18 * BUCKETSIZE + 0)
 #define S_LOWER		(18 * BUCKETSIZE + 1)
@@ -287,10 +265,9 @@ struct tensor {
 #define TESTGT		(19 * BUCKETSIZE + 8)
 #define TESTLE		(19 * BUCKETSIZE + 9)
 #define TESTLT		(19 * BUCKETSIZE + 10)
-#define TGAMMA		(19 * BUCKETSIZE + 11)
-#define TRACE		(19 * BUCKETSIZE + 12)
-#define TRANSPOSE	(19 * BUCKETSIZE + 13)
-#define TTY		(19 * BUCKETSIZE + 14)
+#define TRACE		(19 * BUCKETSIZE + 11)
+#define TRANSPOSE	(19 * BUCKETSIZE + 12)
+#define TTY		(19 * BUCKETSIZE + 13)
 
 #define U_UPPER		(20 * BUCKETSIZE + 0)
 #define U_LOWER		(20 * BUCKETSIZE + 1)
@@ -381,15 +358,19 @@ struct tensor {
 #define BLUE 1
 #define RED 2
 
-#define Trace fprintf(stderr, "%s %d\n", __func__, __LINE__);
+//#define Trace fprintf(stderr, "%s %d\n", __func__, __LINE__);
+#define Trace mp_printf(&mp_plat_print,"[TRACE]%s:%d\n",__func__,__LINE__);
 
-extern struct atom *mem[MAXBLOCKS]; // an array of pointers
+bool noprint = false;
+
+
+extern struct atom *mem; // an array of pointers
 extern struct atom *free_list;
 extern int tos; // top of stack
-extern struct atom *stack[STACKSIZE];
-extern struct atom *symtab[27 * BUCKETSIZE];
-extern struct atom *binding[27 * BUCKETSIZE];
-extern struct atom *usrfunc[27 * BUCKETSIZE];
+extern struct atom **stack;
+extern struct atom **symtab;
+extern struct atom **binding;
+extern struct atom **usrfunc;
 extern struct atom *zero;
 extern struct atom *one;
 extern struct atom *minusone;
@@ -398,11 +379,10 @@ extern int eval_level;
 extern int gc_level;
 extern int expanding;
 extern int drawing;
+extern int nonstop;
 extern int interrupt;
-extern int shuntflag;
-extern int breakflag;
-extern int errorflag;
-extern jmp_buf jmpbuf;
+extern jmp_buf jmpbuf0;
+extern jmp_buf jmpbuf1;
 extern char *trace1;
 extern char *trace2;
 extern int alloc_count;
@@ -421,8 +401,10 @@ extern char strbuf[];
 extern char *outbuf;
 extern int outbuf_index;
 extern int outbuf_length;
+
+
+void init_block(struct atom *mem);
 struct atom * alloc_atom(void);
-void alloc_block(void);
 struct atom * alloc_vector(int nrow);
 struct atom * alloc_matrix(int nrow, int ncol);
 struct atom * alloc_tensor(int nelem);
@@ -464,6 +446,10 @@ int find_denominator(struct atom *p);
 int count_denominators(struct atom *p);
 int count_numerators(struct atom *p);
 void subst(void);
+void evalg(void);
+void evalf(void);
+void evalf_nib(struct atom *p1);
+void evalp(void);
 void eval_abs(struct atom *p1);
 void absfunc(void);
 void eval_add(struct atom *p1);
@@ -502,7 +488,6 @@ void eval_arg(struct atom *p1);
 void argfunc(void);
 void arg_nib(void);
 void eval_binding(struct atom *p1);
-void eval_break(struct atom *p1);
 void eval_ceiling(struct atom *p1);
 void ceilingfunc(void);
 void eval_check(struct atom *p1);
@@ -530,6 +515,7 @@ void dsum(struct atom *p1, struct atom *p2);
 void dproduct(struct atom *p1, struct atom *p2);
 void dpower(struct atom *F, struct atom *X);
 void dlog(struct atom *p1, struct atom *p2);
+void dd(struct atom *p1, struct atom *p2);
 void dfunction(struct atom *p1, struct atom *p2);
 void dsin(struct atom *p1, struct atom *p2);
 void dcos(struct atom *p1, struct atom *p2);
@@ -548,8 +534,6 @@ void derfc(struct atom *p1, struct atom *p2);
 void d_tensor_tensor(struct atom *p1, struct atom *p2);
 void d_scalar_tensor(struct atom *p1, struct atom *p2);
 void d_tensor_scalar(struct atom *p1, struct atom *p2);
-void derivative_of_derivative(struct atom *F, struct atom *X);
-void derivative_of_integral(struct atom *F, struct atom *X);
 void eval_det(struct atom *p1);
 void det(void);
 void eval_dim(struct atom *p1);
@@ -623,26 +607,21 @@ void eval_inner(struct atom *p1);
 void inner(void);
 void eval_integral(struct atom *p1);
 void integral(void);
-void integral_solve(struct atom *F, struct atom *X);
-void integral_solve_nib(struct atom *F, struct atom *X);
+void integral_nib(struct atom *F, struct atom *X);
+void integral_lookup(int h, struct atom *F);
 int integral_classify(struct atom *p);
-int integral_search(int h, struct atom *F, char **table, int n);
+int integral_search(int h, struct atom *F, const char * const *table, int n);
 int integral_search_nib(int h, struct atom *F, struct atom *I, struct atom *C);
 void decomp(void);
 void decomp_sum(struct atom *F, struct atom *X);
 void decomp_product(struct atom *F, struct atom *X);
 void partition_term(void);
-void integral_of_integral(struct atom *F, struct atom *X);
-void integral_of_derivative(struct atom *F, struct atom *X);
 void eval_inv(struct atom *p1);
 void inv(void);
 void eval_kronecker(struct atom *p1);
 void kronecker(void);
-void eval_lgamma(struct atom *p1);
-void lgammafunc(void);
 void eval_log(struct atom *p1);
 void logfunc(void);
-void eval_loop(struct atom *p1);
 void eval_mag(struct atom *p1);
 void magfunc(void);
 void magfunc_nib(void);
@@ -727,7 +706,6 @@ void print_result(void);
 int annotate_result(struct atom *p1, struct atom *p2);
 void eval_product(struct atom *p1);
 void eval_quote(struct atom *p1);
-void eval_rand(struct atom *p1);
 void eval_rank(struct atom *p1);
 void eval_rationalize(struct atom *p1);
 void rationalize(void);
@@ -760,7 +738,7 @@ void set_component(struct atom *LVAL, struct atom *RVAL, int h);
 void setq_usrfunc(struct atom *p1);
 void convert_body(struct atom *A);
 void eval_sgn(struct atom *p1);
-void sgnfunc(void);
+void sgn(void);
 void eval_simplify(struct atom *p1);
 void simplify(void);
 void simplify_nib(void);
@@ -795,18 +773,12 @@ void eval_testgt(struct atom *p1);
 void eval_testle(struct atom *p1);
 void eval_testlt(struct atom *p1);
 int cmp_args(struct atom *p1);
-void eval_tgamma(struct atom *p1);
-void tgammafunc(void);
 void eval_transpose(struct atom *p1);
 void transpose(int n, int m);
 void eval_unit(struct atom *p1);
 void eval_user_function(struct atom *p1);
 void eval_user_symbol(struct atom *p1);
 void eval_zero(struct atom *p1);
-void evalg(void);
-void evalf(void);
-void evalf_nib(struct atom *p1);
-void evalp(void);
 void factor_factor(void);
 void factor_bignum(uint32_t *N, struct atom *M);
 void factor_int(int n);
@@ -857,6 +829,7 @@ void fmt_draw_table(int x, int y, struct atom *p);
 void fmt_putw(uint32_t w);
 void gc(void);
 void untag(struct atom *p);
+//int main(int argc, char *argv[]);
 void run_infile(char *infile);
 void run_stdin(void);
 void display(void);
@@ -898,13 +871,13 @@ int isdenormalpolar(struct atom *p);
 int isdenormalpolarterm(struct atom *p);
 int issquarematrix(struct atom *p);
 int issmallinteger(struct atom *p);
-int dependent(struct atom *f, struct atom *x);
 void run(char *buf);
 void run_buf(char *buf);
 char * scan_input(char *s);
 void print_trace(int color);
 void run_init_script(void);
 void stopf(char *s);
+void exitf(char *s);
 char * scan(char *s);
 char * scan1(char *s);
 char * scan_nib(char *s);
@@ -947,13 +920,60 @@ void set_symbol(struct atom *p1, struct atom *p2, struct atom *p3);
 struct atom * get_binding(struct atom *p1);
 struct atom * get_usrfunc(struct atom *p);
 void init_symbol_table(void);
+
+void eigenmath_init(uint8_t *pHeap,size_t heapSize){
+	eheap_init((void *)pHeap, heapSize);
+	uint32_t sizeOfpAtom = sizeof(struct atom *);
+	uint32_t sizeOfAtom = sizeof(struct atom);
+	//init atom pool
+	MAXATOMS=((heapSize/sizeOfAtom)*4)/10;
+	STACKSIZE=MAXATOMS/12;
+	//mp_printf(&mp_plat_print,"heapsize=%d,mem size = %d\n", heapSize,(uint32_t)(sizeof (struct atom) * MAXATOMS));
+	mem = e_malloc(sizeof (struct atom) * MAXATOMS);
+	if (mem == NULL) {
+		mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("Failed to initialize mem"));
+		return;
+	}
+	init_block(mem);
+	stack = (struct atom **)e_malloc(STACKSIZE * sizeOfpAtom);
+	symtab = (struct atom **)e_malloc((27 * BUCKETSIZE) * sizeOfpAtom);
+	binding = (struct atom **)e_malloc((27 * BUCKETSIZE) * sizeOfpAtom);
+	usrfunc = (struct atom **)e_malloc((27 * BUCKETSIZE) * sizeOfpAtom);
+	zero = NULL;
+	if (stack == NULL || symtab == NULL || binding == NULL || usrfunc == NULL) {
+		mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("Failed to initialize fixed ram area"));
+		return;
+	}
+}
+
+
+
+void init_block(struct atom *mem){
+	//struct atom *p;
+	//p=mem;
+	for (int j = 0; j <MAXATOMS - 1; j++) {
+		mem[j].atomtype = FREEATOM;
+		mem[j].u.next = &mem [j + 1];
+	}
+	mem[MAXATOMS - 1].atomtype = FREEATOM;
+	mem[MAXATOMS - 1].u.next = NULL;
+	free_list = mem;
+	free_count = MAXATOMS;
+	alloc_count = 0;
+ // return the first block pointer
+}
+
 struct atom *
 alloc_atom(void)
 {
 	struct atom *p;
 
-	if (free_count == 0)
-		alloc_block();
+	if (free_count == 0){
+		tos = 0;
+		gc(); // prep for next run
+		exitf("out of atoms pool");
+	}
+		
 
 	p = free_list;
 	free_list = p->u.next;
@@ -964,34 +984,6 @@ alloc_atom(void)
 	return p;
 }
 
-void
-alloc_block(void)
-{
-	int i;
-	struct atom *p;
-
-	if (block_count == MAXBLOCKS) {
-		tos = 0;
-		gc(); // prep for next run
-		stopf("out of memory");
-	}
-
-	p = alloc_mem(BLOCKSIZE * sizeof (struct atom));
-
-	mem[block_count++] = p;
-
-	for (i = 0; i < BLOCKSIZE - 1; i++) {
-		p[i].atomtype = FREEATOM;
-		p[i].u.next = p + i + 1;
-	}
-
-	p[i].atomtype = FREEATOM;
-	p[i].u.next = NULL;
-
-	free_list = p;
-	free_count = BLOCKSIZE;
-}
-
 struct atom *
 alloc_vector(int nrow)
 {
@@ -1042,9 +1034,9 @@ alloc_str(void)
 void *
 alloc_mem(int n)
 {
-	void *p = m_malloc(n);
+	void *p = e_malloc(n);
 	if (p == NULL)
-		exit(1);
+		stopf("alloc mem failed out of memory");
 	return p;
 }
 double
@@ -1079,7 +1071,7 @@ mscan(char *s)
 	int i, k, len;
 	uint32_t *a, *b, *t;
 	a = mint(0);
-	len = (int) strlen(s);
+	len = (int) strlen((const char *)s);
 	if (len == 0)
 		return a;
 	k = len % 9;
@@ -1127,7 +1119,7 @@ mstr(uint32_t *u)
 
 	if (n > len) {
 		if (buf)
-			m_free(buf);
+			e_free(buf);
 		buf = alloc_mem(n);
 		len = n;
 	}
@@ -1509,7 +1501,7 @@ mnew(int n)
 void
 mfree(uint32_t *u)
 {
-	m_free(u - 1);
+	e_free(u - 1);
 	bignum_count--;
 }
 
@@ -1979,7 +1971,7 @@ subst(void)
 void
 evalg(void)
 {
-	if (gc_level == eval_level && alloc_count > MAXBLOCKS * BLOCKSIZE / 10)
+	if (gc_level == eval_level && alloc_count > MAXATOMS / 10)
 		gc();
 	gc_level++;
 	evalf();
@@ -2008,10 +2000,8 @@ evalf_nib(struct atom *p1)
 	if (interrupt)
 		stopf("interrupt");
 
-	// this is needed to prevent seg fault (STACKSIZE is greater than process stack)
-
-	if (eval_level > 1000)
-		stopf("evaluation depth exceeded, possibly due to recursive function or circular symbol definition");
+	if (eval_level == 200)
+		stopf("circular definition?");
 
 	if (eval_level > max_eval_level)
 		max_eval_level = eval_level;
@@ -3512,12 +3502,6 @@ eval_binding(struct atom *p1)
 	push(get_binding(cadr(p1)));
 }
 void
-eval_break(struct atom *p1)
-{
-	breakflag = 1;
-	push_symbol(NIL);
-}
-void
 eval_ceiling(struct atom *p1)
 {
 	push(cadr(p1));
@@ -4202,9 +4186,6 @@ eval_defint(struct atom *p1)
 		integral();
 		F = pop();
 
-		if (findf(F, symbol(INTEGRAL)))
-			stopf("defint: unsolved integral");
-
 		push(F);
 		push(X);
 		push(B);
@@ -4333,16 +4314,6 @@ d_scalar_scalar(struct atom *F, struct atom *X)
 	if (!isusersymbol(X))
 		stopf("derivative: symbol expected");
 
-	if (car(F) == symbol(DERIVATIVE)) {
-		derivative_of_derivative(F, X);
-		return;
-	}
-
-	if (car(F) == symbol(INTEGRAL)) {
-		derivative_of_integral(F, X);
-		return;
-	}
-
 	// d(x,x)?
 
 	if (equal(F, X)) {
@@ -4372,6 +4343,11 @@ d_scalar_scalar(struct atom *F, struct atom *X)
 		return;
 	}
 
+	if (car(F) == symbol(DERIVATIVE)) {
+		dd(F, X);
+		return;
+	}
+
 	if (car(F) == symbol(LOG)) {
 		dlog(F, X);
 		return;
@@ -4447,6 +4423,11 @@ d_scalar_scalar(struct atom *F, struct atom *X)
 		return;
 	}
 
+	if (car(F) == symbol(INTEGRAL) && caddr(F) == X) {
+		push(cadr(F));
+		return;
+	}
+
 	dfunction(F, X);
 }
 
@@ -4542,6 +4523,58 @@ dlog(struct atom *p1, struct atom *p2)
 	divide();
 }
 
+//	derivative of derivative
+//
+//	example: d(d(f(x,y),y),x)
+//
+//	p1 = d(f(x,y),y)
+//
+//	p2 = x
+//
+//	cadr(p1) = f(x,y)
+//
+//	caddr(p1) = y
+
+void
+dd(struct atom *p1, struct atom *p2)
+{
+	struct atom *p3;
+
+	// d(f(x,y),x)
+
+	push(cadr(p1));
+	push(p2);
+	derivative();
+
+	p3 = pop();
+
+	if (car(p3) == symbol(DERIVATIVE)) {
+
+		// sort dx terms
+
+		push_symbol(DERIVATIVE);
+		push_symbol(DERIVATIVE);
+		push(cadr(p3));
+
+		if (lessp(caddr(p3), caddr(p1))) {
+			push(caddr(p3));
+			list(3);
+			push(caddr(p1));
+		} else {
+			push(caddr(p1));
+			list(3);
+			push(caddr(p3));
+		}
+
+		list(3);
+
+	} else {
+		push(p3);
+		push(caddr(p1));
+		derivative();
+	}
+}
+
 // derivative of a generic function
 
 void
@@ -4855,79 +4888,6 @@ d_tensor_scalar(struct atom *p1, struct atom *p2)
 
 	push(p3);
 }
-
-void
-derivative_of_derivative(struct atom *F, struct atom *X)
-{
-	struct atom *G, *Y;
-
-	G = cadr(F);
-	Y = caddr(F);
-
-	if (X == Y) {
-		push_symbol(DERIVATIVE);
-		push(F);
-		push(X);
-		list(3);
-		return;
-	}
-
-	push(G);
-	push(X);
-	derivative();
-
-	G = pop();
-
-	if (car(G) == symbol(DERIVATIVE)) {
-
-		// sort derivatives
-
-		F = cadr(G);
-		X = caddr(G);
-
-		push_symbol(DERIVATIVE);
-		push_symbol(DERIVATIVE);
-		push(F);
-
-		if (lessp(X, Y)) {
-			push(X);
-			list(3);
-			push(Y);
-			list(3);
-		} else {
-			push(Y);
-			list(3);
-			push(X);
-			list(3);
-		}
-
-		return;
-	}
-
-	push(G);
-	push(Y);
-	derivative();
-}
-
-void
-derivative_of_integral(struct atom *F, struct atom *X)
-{
-	struct atom *G, *Y;
-
-	G = cadr(F);
-	Y = caddr(F);
-
-	if (X == Y) {
-		push(G); // derivative and integral cancel
-		return;
-	}
-
-	push(G);
-	push(X);
-	derivative();
-	push(Y);
-	integral();
-}
 void
 eval_det(struct atom *p1)
 {
@@ -5095,9 +5055,9 @@ eval_eigenvec(struct atom *p1)
 				stopf("eigenvec");
 
 	if (D)
-		m_free(D);
+		e_free(D);
 	if (Q)
-		m_free(Q);
+		e_free(Q);
 
 	D = alloc_mem(n * n * sizeof (double));
 	Q = alloc_mem(n * n * sizeof (double));
@@ -5420,8 +5380,6 @@ eigenvec_step_nib(double *D, double *Q, int n, int p, int q)
 	D[n * p + q] = 0.0;
 	D[n * q + p] = 0.0;
 }
-// first author: github.com/phbillet
-
 void
 eval_erf(struct atom *p1)
 {
@@ -5435,7 +5393,7 @@ erffunc(void)
 {
 	int i, n;
 	double d;
-	struct atom *p1, *p2;
+	struct atom *p1;
 
 	p1 = pop();
 
@@ -5451,32 +5409,32 @@ erffunc(void)
 		return;
 	}
 
-	push(p1);
-	floatfunc();
-	p2 = pop();
-
-	if (isnum(p2)) {
-		push(p2);
+	if (isdouble(p1)) {
+		push(p1);
 		d = pop_double();
 		d = erf(d);
 		push_double(d);
 		return;
 	}
 
+	if (iszero(p1)) {
+		push_integer(0);
+		return;
+	}
+
 	if (isnegativeterm(p1)) {
 		push_symbol(ERF);
 		push(p1);
 		negate();
 		list(2);
 		negate();
-	} else {
-		push_symbol(ERF);
-		push(p1);
-		list(2);
+		return;
 	}
-}
-// first author: github.com/phbillet
 
+	push_symbol(ERF);
+	push(p1);
+	list(2);
+}
 void
 eval_erfc(struct atom *p1)
 {
@@ -5490,7 +5448,7 @@ erfcfunc(void)
 {
 	int i, n;
 	double d;
-	struct atom *p1, *p2;
+	struct atom *p1;
 
 	p1 = pop();
 
@@ -5506,32 +5464,22 @@ erfcfunc(void)
 		return;
 	}
 
-	push(p1);
-	floatfunc();
-	p2 = pop();
-
-	if (isnum(p2)) {
-		push(p2);
+	if (isdouble(p1)) {
+		push(p1);
 		d = pop_double();
 		d = erfc(d);
 		push_double(d);
 		return;
 	}
 
-	if (isnegativeterm(p1)) {
-		push_symbol(ERF);
-		push(p1);
-		negate();
-		list(2);
-	} else {
-		push_symbol(ERF);
-		push(p1);
-		list(2);
-		negate();
+	if (iszero(p1)) {
+		push_integer(1);
+		return;
 	}
 
-	push_integer(1);
-	add();
+	push_symbol(ERFC);
+	push(p1);
+	list(2);
 }
 void
 eval_eval(struct atom *p1)
@@ -6225,7 +6173,7 @@ floorfunc(void)
 void
 eval_for(struct atom *p1)
 {
-	int j, k, t;
+	int j, k;
 	struct atom *p2, *p3;
 
 	p2 = cadr(p1);
@@ -6252,9 +6200,6 @@ eval_for(struct atom *p1)
 
 	save_symbol(p2);
 
-	t = breakflag;
-	breakflag = 0;
-
 	for (;;) {
 		push_integer(j);
 		p3 = pop();
@@ -6263,14 +6208,8 @@ eval_for(struct atom *p1)
 		while (iscons(p3)) {
 			push(car(p3));
 			evalg();
-			pop();
+			pop(); // discard return value
 			p3 = cdr(p3);
-			if (breakflag || errorflag) {
-				breakflag = t;
-				restore_symbol();
-				push_symbol(NIL);
-				return;
-			}
 		}
 		if (j == k)
 			break;
@@ -6280,9 +6219,9 @@ eval_for(struct atom *p1)
 			j--;
 	}
 
-	breakflag = t;
 	restore_symbol();
-	push_symbol(NIL);
+
+	push_symbol(NIL); // return value
 }
 void
 eval_hadamard(struct atom *p1)
@@ -6439,13 +6378,8 @@ indexfunc(struct atom *T, int h)
 	for (i = 0; i < n; i++) {
 		push(stack[h + i]);
 		t = pop_integer();
-		if (t < 1 || t > T->u.tensor->dim[i]) {
-			if (shuntflag) {
-				errorflag = 1;
-				t = 1;
-			} else
-				stopf("index error");
-		}
+		if (t < 1 || t > T->u.tensor->dim[i])
+			stopf("index error");
 		k = k * T->u.tensor->dim[i] + t - 1;
 	}
 
@@ -6496,7 +6430,10 @@ print_infixform(struct atom *p)
 	outbuf_init();
 	infixform_expr(p);
 	outbuf_puts("\n");
-	printbuf(outbuf, BLACK);
+	outbuf_puts("\0");
+	if (noprint == false){
+		printbuf(outbuf, BLACK);
+	}	
 }
 
 void
@@ -7127,7 +7064,7 @@ inner(void)
 
 	push(p3);
 }
-char *integral_tab_exp[] = {
+const char * const integral_tab_exp[] = {
 
 // x^n exp(a x + b)
 
@@ -7370,7 +7307,7 @@ char *integral_tab_exp[] = {
 
 // log(a x) is transformed to log(a) + log(x)
 
-char *integral_tab_log[] = {
+const char * const integral_tab_log[] = {
 
 	"log(x)",
 	"x log(x) - x",
@@ -7417,7 +7354,7 @@ char *integral_tab_log[] = {
 	"1",
 };
 
-char *integral_tab_trig[] = {
+const char * const integral_tab_trig[] = {
 
 	"sin(a x)",
 	"-cos(a x) / a",
@@ -7676,7 +7613,7 @@ char *integral_tab_trig[] = {
 	"1",
 };
 
-char *integral_tab_power[] = {
+const char * const integral_tab_power[] = {
 
 	"a", // for forms c^d where both c and d are constant expressions
 	"a x",
@@ -7797,7 +7734,7 @@ char *integral_tab_power[] = {
 	"1",
 };
 
-char *integral_tab[] = {
+const char * const integral_tab[] = {
 
 	"a",
 	"a x",
@@ -8162,34 +8099,37 @@ integral(void)
 		push(X);
 		partition_term();	// push const part then push var part
 		F = pop();		// pop var part
-		integral_solve(F, X);
+		integral_nib(F, X);
 		multiply();		// multiply by const part
 		return;
 	}
 
-	integral_solve(F, X);
+	integral_nib(F, X);
 }
 
 void
-integral_solve(struct atom *F, struct atom *X)
+integral_nib(struct atom *F, struct atom *X)
 {
+	int h;
 	struct atom *p;
 
-	if (car(F) == symbol(INTEGRAL)) {
-		integral_of_integral(F, X);
-		return;
-	}
-
-	if (car(F) == symbol(DERIVATIVE)) {
-		integral_of_derivative(F, X);
-		return;
-	}
-
 	save_symbol(symbol(SA));
 	save_symbol(symbol(SB));
 	save_symbol(symbol(SX));
 
-	integral_solve_nib(F, X);
+	set_symbol(symbol(SX), X, symbol(NIL));
+
+	// put constants in F(X) on the stack
+
+	h = tos;
+
+	push_integer(1); // 1 is a candidate for a or b
+
+	push(F);
+	push(X);
+	decomp(); // push const coeffs
+
+	integral_lookup(h, F);
 
 	p = pop();
 	restore_symbol();
@@ -8199,19 +8139,9 @@ integral_solve(struct atom *F, struct atom *X)
 }
 
 void
-integral_solve_nib(struct atom *F, struct atom *X)
+integral_lookup(int h, struct atom *F)
 {
-	int h, t;
-
-	set_symbol(symbol(SX), X, symbol(NIL));
-
-	h = tos;
-
-	push_integer(1); // 1 is a candidate for a and b
-
-	push(F);
-	push(X);
-	decomp(); // push possible substitutions for a and b
+	int t;
 
 	t = integral_classify(F);
 
@@ -8232,12 +8162,7 @@ integral_solve_nib(struct atom *F, struct atom *X)
 			return;
 	}
 
-	tos = h; // pop all
-
-	push_symbol(INTEGRAL);
-	push(F);
-	push(X);
-	list(3);
+	stopf("integral: no solution found");
 }
 
 int
@@ -8266,17 +8191,17 @@ integral_classify(struct atom *p)
 }
 
 int
-integral_search(int h, struct atom *F, char **table, int n)
+integral_search(int h, struct atom *F, const char * const *table, int n)
 {
 	int i;
 	struct atom *C, *I;
 
 	for (i = 0; i < n; i += 3) {
 
-		scan1(table[i + 0]); // integrand
+		scan1((char *)table[i + 0]); // integrand
 		I = pop();
 
-		scan1(table[i + 2]); // condition
+		scan1((char *)table[i + 2]); // condition
 		C = pop();
 
 		if (integral_search_nib(h, F, I, C))
@@ -8288,7 +8213,7 @@ integral_search(int h, struct atom *F, char **table, int n)
 
 	tos = h; // pop all
 
-	scan1(table[i + 1]); // answer
+	scan1((char *)table[i + 1]); // answer
 	evalf();
 
 	return 1;
@@ -8339,7 +8264,7 @@ decomp(void)
 
 	// is the entire expression constant?
 
-	if (!dependent(F, X)) {
+	if (!findf(F, X)) {
 		push(F);
 		return;
 	}
@@ -8383,7 +8308,7 @@ decomp_sum(struct atom *F, struct atom *X)
 
 	while (iscons(p1)) {
 		p2 = car(p1);
-		if (dependent(p2, X)) {
+		if (findf(p2, X)) {
 			if (car(p2) == symbol(MULTIPLY)) {
 				push(p2);
 				push(X);
@@ -8433,7 +8358,7 @@ decomp_sum(struct atom *F, struct atom *X)
 	h = tos;
 	p1 = cdr(F);
 	while (iscons(p1)) {
-		if (!dependent(car(p1), X))
+		if (!findf(car(p1), X))
 			push(car(p1));
 		p1 = cdr(p1);
 	}
@@ -8458,7 +8383,7 @@ decomp_product(struct atom *F, struct atom *X)
 
 	p1 = cdr(F);
 	while (iscons(p1)) {
-		if (dependent(car(p1), X)) {
+		if (findf(car(p1), X)) {
 			push(car(p1));
 			push(X);
 			decomp();
@@ -8471,7 +8396,7 @@ decomp_product(struct atom *F, struct atom *X)
 	h = tos;
 	p1 = cdr(F);
 	while (iscons(p1)) {
-		if (!dependent(car(p1), X))
+		if (!findf(car(p1), X))
 			push(car(p1));
 		p1 = cdr(p1);
 	}
@@ -8500,7 +8425,7 @@ partition_term(void)
 	h = tos;
 	p1 = cdr(F);
 	while (iscons(p1)) {
-		if (!dependent(car(p1), X))
+		if (!findf(car(p1), X))
 			push(car(p1));
 		p1 = cdr(p1);
 	}
@@ -8521,7 +8446,7 @@ partition_term(void)
 	h = tos;
 	p1 = cdr(F);
 	while (iscons(p1)) {
-		if (dependent(car(p1), X))
+		if (findf(car(p1), X))
 			push(car(p1));
 		p1 = cdr(p1);
 	}
@@ -8537,100 +8462,6 @@ partition_term(void)
 		cons(); // makes MULTIPLY head of list
 	}
 }
-
-void
-integral_of_integral(struct atom *F, struct atom *X)
-{
-	struct atom *G, *Y;
-
-	G = cadr(F);
-	Y = caddr(F);
-
-	// if X == Y then F is not integrable for X
-
-	if (equal(X, Y)) {
-		push_symbol(INTEGRAL);
-		push(F);
-		push(X);
-		list(3);
-		return;
-	}
-
-	push(G);
-	push(X);
-	integral();
-
-	G = pop();
-
-	if (car(G) == symbol(INTEGRAL)) {
-
-		// sort integrals by measure
-
-		F = cadr(G);
-		X = caddr(G);
-
-		push_symbol(INTEGRAL);
-		push_symbol(INTEGRAL);
-		push(F);
-
-		if (lessp(X, Y)) {
-			push(X);
-			list(3);
-			push(Y);
-			list(3);
-		} else {
-			push(Y);
-			list(3);
-			push(X);
-			list(3);
-		}
-
-		return;
-	}
-
-	push(G);
-	push(Y);
-	integral();
-}
-
-void
-integral_of_derivative(struct atom *F, struct atom *X)
-{
-	struct atom *G, *Y;
-
-	G = cadr(F);
-	Y = caddr(F);
-
-	if (equal(X, Y)) {
-		push(G); // integral and derivative cancel
-		return;
-	}
-
-	push(G);
-	push(X);
-	integral();
-
-	G = pop();
-
-	// integral before derivative
-
-	if (car(G) == symbol(INTEGRAL)) {
-		F = cadr(G);
-		X = caddr(G);
-		push_symbol(INTEGRAL);
-		push_symbol(DERIVATIVE);
-		push(F);
-		push(Y);
-		list(3);
-		push(X);
-		list(3);
-		return;
-	}
-
-	push(G);
-	push(Y);
-	derivative();
-}
 void
 eval_inv(struct atom *p1)
 {
@@ -8730,51 +8561,6 @@ kronecker(void)
 	push(p3);
 }
 void
-eval_lgamma(struct atom *p1)
-{
-	push(cadr(p1));
-	evalf();
-	lgammafunc();
-}
-
-void
-lgammafunc(void)
-{
-	int i, n;
-	double d;
-	struct atom *p1, *p2;
-
-	p1 = pop();
-
-	if (istensor(p1)) {
-		p1 = copy_tensor(p1);
-		n = p1->u.tensor->nelem;
-		for (i = 0; i < n; i++) {
-			push(p1->u.tensor->elem[i]);
-			lgammafunc();
-			p1->u.tensor->elem[i] = pop();
-		}
-		push(p1);
-		return;
-	}
-
-	push(p1);
-	floatfunc();
-	p2 = pop();
-
-	if (isnum(p2)) {
-		push(p2);
-		d = pop_double();
-		d = lgamma(d);
-		push_double(d);
-		return;
-	}
-
-	push_symbol(LGAMMA);
-	push(p1);
-	list(2);
-}
-void
 eval_log(struct atom *p1)
 {
 	push(cadr(p1));
@@ -8918,32 +8704,6 @@ logfunc(void)
 	list(2);
 }
 void
-eval_loop(struct atom *p1)
-{
-	int t;
-	struct atom *p2;
-	t = breakflag;
-	breakflag = 0;
-	if (lengthf(p1) < 2) {
-		push_symbol(NIL);
-		return;
-	}
-	for (;;) {
-		p2 = cdr(p1);
-		while (iscons(p2)) {
-			push(car(p2));
-			evalg();
-			pop();
-			p2 = cdr(p2);
-			if (breakflag || errorflag) {
-				breakflag = t;
-				push_symbol(NIL);
-				return;
-			}
-		}
-	}
-}
-void
 eval_mag(struct atom *p1)
 {
 	push(cadr(p1));
@@ -9533,9 +9293,8 @@ void
 expand_sum_factors(int h)
 {
 	int i, n;
-	struct atom *p1, *p2;
-	// Here, compiler treats warning as error, we use NULL to supress compiler error.
-	p2 = NULL;
+	struct atom *p1=NULL;
+	struct atom *p2=NULL;
 
 	if (tos - h < 2)
 		return;
@@ -9962,9 +9721,9 @@ nroots(void)
 	n = tos - h; // number of coeffs on stack
 
 	if (cr)
-		m_free(cr);
+		e_free(cr);
 	if (ci)
-		m_free(ci);
+		e_free(ci);
 
 	cr = alloc_mem(n * sizeof (double));
 	ci = alloc_mem(n * sizeof (double));
@@ -10627,12 +10386,8 @@ power_numbers(struct atom *BASE, struct atom *EXPO)
 	// 0^n
 
 	if (iszero(BASE)) {
-		if (isnegativenumber(EXPO)) {
-			if (shuntflag)
-				errorflag = 1;
-			else
-				stopf("divide by zero");
-		}
+		if (isnegativenumber(EXPO))
+			stopf("divide by zero");
 		push_integer(0);
 		return;
 	}
@@ -11680,7 +11435,10 @@ print_prefixform(struct atom *p)
 	outbuf_init();
 	prefixform(p);
 	outbuf_puts("\n");
-	printbuf(outbuf, BLACK);
+	outbuf_puts("\0");
+	if (noprint == false){
+		printbuf(outbuf, BLACK);
+	}
 }
 
 void
@@ -11722,8 +11480,9 @@ prefixform(struct atom *p)
 	case DOUBLE:
 		snprintf(strbuf, STRBUFLEN, "%g", p->u.d);
 		outbuf_puts(strbuf);
-		if (!strchr(strbuf, '.') && !strchr(strbuf, 'e'))
+		if (!strchr(strbuf, '.') && !strchr(strbuf, 'e')){
 			outbuf_puts(".0");
+		}
 		break;
 	case KSYM:
 	case USYM:
@@ -11775,8 +11534,10 @@ print_result(void)
 	if (p1 == symbol(TTY) || iszero(p1)) {
 		push(p2);
 		display();
-	} else
+	} else{
 		print_infixform(p2);
+	}
+		
 }
 
 // returns 1 if result should be annotated
@@ -11875,13 +11636,6 @@ eval_quote(struct atom *p1)
 	push(cadr(p1)); // not evaluated
 }
 void
-eval_rand(struct atom *p1)
-{
-	double d;
-	d = (double) rand() / ((double) RAND_MAX + 1);
-	push_double(d);
-}
-void
 eval_rank(struct atom *p1)
 {
 	push(cadr(p1));
@@ -12738,7 +12492,7 @@ read_file(char *filename)
 	}
 
 	if (read(fd, buf, n) != n) {
-		m_free(buf);
+		e_free(buf);
 		close(fd);
 		return NULL;
 	}
@@ -12994,21 +12748,19 @@ convert_body(struct atom *A)
 	push_symbol(ARG9);
 	subst();
 }
-// first author: github.com/phbillet
-
 void
 eval_sgn(struct atom *p1)
 {
 	push(cadr(p1));
 	evalf();
-	sgnfunc();
+	sgn();
 }
 
 void
-sgnfunc(void)
+sgn(void)
 {
 	int i, n;
-	struct atom *p1, *p2;
+	struct atom *p1;
 
 	p1 = pop();
 
@@ -13017,34 +12769,26 @@ sgnfunc(void)
 		n = p1->u.tensor->nelem;
 		for (i = 0; i < n; i++) {
 			push(p1->u.tensor->elem[i]);
-			sgnfunc();
+			sgn();
 			p1->u.tensor->elem[i] = pop();
 		}
 		push(p1);
 		return;
 	}
 
-	p2 = p1;
-
-	if (!isnum(p2)) {
-		push(p2);
-		floatfunc();
-		p2 = pop();
-	}
-
-	if (!isnum(p2)) {
+	if (!isnum(p1)) {
 		push_symbol(SGN);
 		push(p1);
 		list(2);
 		return;
 	}
 
-	if (iszero(p2)) {
+	if (iszero(p1)) {
 		push_integer(0);
 		return;
 	}
 
-	if (isnegativenumber(p2))
+	if (isnegativenumber(p1))
 		push_integer(-1);
 	else
 		push_integer(1);
@@ -13623,8 +13367,6 @@ eval_status(struct atom *p1)
 
 	outbuf_init();
 
-	snprintf(strbuf, STRBUFLEN, "block_count %d (%d%%)\n", block_count, 100 * block_count / MAXBLOCKS);
-	outbuf_puts(strbuf);
 
 	snprintf(strbuf, STRBUFLEN, "free_count %d\n", free_count);
 	outbuf_puts(strbuf);
@@ -13650,10 +13392,11 @@ eval_status(struct atom *p1)
 	snprintf(strbuf, STRBUFLEN, "max_eval_level %d\n", max_eval_level);
 	outbuf_puts(strbuf);
 
-	snprintf(strbuf, STRBUFLEN, "max_tos %d (%d%%)\n", max_tos, 100 * max_tos / STACKSIZE);
+	snprintf(strbuf, STRBUFLEN, "max_tos %d (%ld%%)\n", max_tos, (int32_t)(100 * max_tos / STACKSIZE));
 	outbuf_puts(strbuf);
-
-	printbuf(outbuf, BLACK);
+	if (noprint == false){
+		printbuf(outbuf, BLACK);
+	}
 
 	push_symbol(NIL);
 }
@@ -14276,58 +14019,13 @@ cmp_args(struct atom *p1)
 	if (iszero(p1))
 		return 0;
 	if (!isnum(p1))
-		stopf("arithmetic comparison: not a number");
+		stopf("compare err");
 	if (isnegativenumber(p1))
 		return -1;
 	else
 		return 1;
 }
 void
-eval_tgamma(struct atom *p1)
-{
-	push(cadr(p1));
-	evalf();
-	tgammafunc();
-}
-
-void
-tgammafunc(void)
-{
-	int i, n;
-	double d;
-	struct atom *p1, *p2;
-
-	p1 = pop();
-
-	if (istensor(p1)) {
-		p1 = copy_tensor(p1);
-		n = p1->u.tensor->nelem;
-		for (i = 0; i < n; i++) {
-			push(p1->u.tensor->elem[i]);
-			tgammafunc();
-			p1->u.tensor->elem[i] = pop();
-		}
-		push(p1);
-		return;
-	}
-
-	push(p1);
-	floatfunc();
-	p2 = pop();
-
-	if (isnum(p2)) {
-		push(p2);
-		d = pop_double();
-		d = tgamma(d);
-		push_double(d);
-		return;
-	}
-
-	push_symbol(TGAMMA);
-	push(p1);
-	list(2);
-}
-void
 eval_transpose(struct atom *p1)
 {
 	int m, n;
@@ -14471,7 +14169,7 @@ eval_user_function(struct atom *p1)
 		push(FUNC_NAME);
 		while (iscons(FUNC_ARGS)) {
 			push(car(FUNC_ARGS));
-			evalg();
+			evalg(); // p1 is on frame stack, not reclaimed
 			FUNC_ARGS = cdr(FUNC_ARGS);
 		}
 		list(tos - h);
@@ -14701,7 +14399,7 @@ factor_bignum(uint32_t *N, struct atom *M)
 
 #define NPRIME 4792
 
-int primetab[NPRIME] = {
+const int primetab[NPRIME] = {
 2,3,5,7,11,13,17,19,
 23,29,31,37,41,43,47,53,
 59,61,67,71,73,79,83,89,
@@ -15321,7 +15019,7 @@ factor_int(int n)
 
 	for (k = 0; k < NPRIME; k++) {
 
-		d = primetab[k];
+		d = (int)primetab[k];
 
 		m = 0;
 
@@ -15364,11 +15062,11 @@ factor_int(int n)
 #define VAL2(p) ((int) cadr(p)->u.d)
 
 #define PLUS_SIGN '+'
-#define MINUS_SIGN 0xe28892
-#define MULTIPLY_SIGN 0xc397
-#define GREATEREQUAL 0xe289a5
-#define LESSEQUAL 0xe289a4
-
+#define MINUS_SIGN 0x2D       // '-' 0xe28892
+#define MULTIPLY_SIGN 0x2A       // '*' 0xc397
+#define GREATEREQUAL 0xF2 //   IBM code page 437   0xe289a5
+#define LESSEQUAL 0xF3//   IBM code page 437   0xe289a4
+/*
 #define BDLL 0xe295b4 // BOX DRAW LIGHT LEFT
 #define BDLR 0xe295b6 // BOX DRAW LIGHT RIGHT
 
@@ -15379,9 +15077,23 @@ factor_int(int n)
 #define BDLDAL 0xe29490 // BOX DRAW LIGHT DOWN AND LEFT
 #define BDLUAR 0xe29494 // BOX DRAW LIGHT UP AND RIGHT
 #define BDLUAL 0xe29498 // BOX DRAW LIGHT UP AND LEFT
+*/
+/*  ─ */
+#define BDLL      0xC4   /* ╴→ ─ */
+#define BDLR      0xC4   /* ╶→ ─ */
 
-// #define MAX(a,b) ((a) > (b) ? (a) : (b))
-// #define MIN(a,b) ((a) < (b) ? (a) : (b))
+/*  */
+#define BDLH      0xC4   /* ─ */
+#define BDLV      0xB3   /* │ */
+
+/*  */
+#define BDLDAR    0xDA   /* ┌ */
+#define BDLDAL    0xBF   /* ┐ */
+#define BDLUAR    0xC0   /* └ */
+#define BDLUAL    0xD9   /* ┘ */
+
+//#define MAX(a,b) ((a) > (b) ? (a) : (b))
+//#define MIN(a,b) ((a) < (b) ? (a) : (b))
 
 int fmt_level;
 int fmt_nrow;
@@ -15416,8 +15128,11 @@ fmt(void)
 
 	if (m > fmt_buf_len) {
 		if (fmt_buf)
-			m_free(fmt_buf);
+			e_free(fmt_buf);
 		fmt_buf = alloc_mem(m);
+		if (!fmt_buf) {
+			stopf("fmt: out of memory allocating %d bytes");
+		}
 		fmt_buf_len = m;
 	}
 
@@ -15436,8 +15151,11 @@ fmt(void)
 	}
 
 	fmt_putw('\n'); // blank line after result
-
-	printbuf(outbuf, BLACK);
+	//fmt_putw('\0'); 
+	if (noprint==false){
+		printbuf(outbuf, BLACK);
+	}
+	
 }
 
 void
@@ -15735,11 +15453,11 @@ fmt_function(struct atom *p)
 
 	// default
 
-	if (issymbol(car(p)))
+	if (issymbol(car(p))){
 		fmt_symbol(car(p));
-	else
+	}else{
 		fmt_subexpr(car(p));
-
+	}
 	fmt_args(p);
 }
 
@@ -15793,16 +15511,18 @@ fmt_matrix(struct atom *p, int d, int k)
 
 	span = 1;
 
-	for (i = d + 2; i < p->u.tensor->ndim; i++)
+	for (i = d + 2; i < p->u.tensor->ndim; i++){
 		span *= p->u.tensor->dim[i];
+	}
 
 	n = p->u.tensor->dim[d];	// number of rows
 	m = p->u.tensor->dim[d + 1];	// number of columns
 
-	for (i = 0; i < n; i++)
-		for (j = 0; j < m; j++)
+	for (i = 0; i < n; i++){
+		for (j = 0; j < m; j++){
 			fmt_matrix(p, d + 2, k + (i * m + j) * span);
-
+		}
+	}
 	fmt_update_table(n, m);
 }
 
@@ -15822,11 +15542,13 @@ fmt_numerators(struct atom *p)
 		q = car(p);
 		p = cdr(p);
 
-		if (!isnumerator(q))
+		if (!isnumerator(q)){
 			continue;
+		}
 
-		if (tos > t)
+		if (tos > t){
 			fmt_space();
+		}
 
 		if (isrational(q)) {
 			s = mstr(q->u.q.a);
@@ -15834,15 +15556,16 @@ fmt_numerators(struct atom *p)
 			continue;
 		}
 
-		if (car(q) == symbol(ADD) && n == 1)
+		if (car(q) == symbol(ADD) && n == 1){
 			fmt_expr(q); // parens not needed
-		else
+		}else{
 			fmt_factor(q);
+		}
 	}
 
-	if (t == tos)
+	if (t == tos){
 		fmt_roman_char('1'); // no numerators
-
+	}
 	fmt_update_list(t);
 }
 
@@ -15945,9 +15668,9 @@ fmt_reciprocal(struct atom *p)
 
 	t = tos;
 
-	if (isminusone(caddr(p)))
+	if (isminusone(caddr(p))){
 		fmt_expr(cadr(p));
-	else {
+	}else {
 		fmt_base(cadr(p));
 		fmt_numeric_exponent(caddr(p)); // sign is not emitted
 	}
@@ -15978,8 +15701,9 @@ fmt_roman_char(int c)
 void
 fmt_roman_string(char *s)
 {
-	while (*s)
+	while (*s){
 		fmt_roman_char(*s++);
+	}
 }
 
 void
@@ -16026,8 +15750,9 @@ fmt_symbol(struct atom *p)
 
 	k = fmt_symbol_fragment(s, 0);
 
-	if (s[k] == '\0')
+	if (s[k] == '\0'){
 		return;
+	}
 
 	// emit subscript
 
@@ -16035,8 +15760,9 @@ fmt_symbol(struct atom *p)
 
 	t = tos;
 
-	while (s[k] != '\0')
+	while (s[k] != '\0'){
 		k = fmt_symbol_fragment(s, k);
+	}
 
 	fmt_update_list(t);
 
@@ -16047,7 +15773,7 @@ fmt_symbol(struct atom *p)
 
 #define NUM_SYMBOL_NAMES 49
 
-char *symbol_name_tab[NUM_SYMBOL_NAMES] = {
+const char * const symbol_name_tab[NUM_SYMBOL_NAMES] = {
 
 	"Alpha",
 	"Beta",
@@ -16102,57 +15828,57 @@ char *symbol_name_tab[NUM_SYMBOL_NAMES] = {
 	"hbar",
 };
 
-int symbol_unicode_tab[NUM_SYMBOL_NAMES] = {
+const int symbol_unicode_tab[NUM_SYMBOL_NAMES] = {
 
-	0xce91, // Alpha
-	0xce92, // Beta
-	0xce93, // Gamma
-	0xce94, // Delta
-	0xce95, // Epsilon
-	0xce96, // Zeta
-	0xce97, // Eta
-	0xce98, // Theta
-	0xce99, // Iota
-	0xce9a, // Kappa
+	0x41, // Alpha  0xce91
+	0x42, // Beta 0xce92
+	0xE2, // Gamma 0xce93
+	0x7F, // Delta0xce94
+	0x45, // Epsilon0xce95
+	0x5A, // Zeta0xce96
+	0x48, // Eta0xce97
+	0xE9, // Theta0xce98
+	0x49, // Iota0xce99
+	0x4B, // Kappa0xce9a
 	0xce9b, // Lambda
-	0xce9c, // Mu
-	0xce9d, // Nu
+	0x4D, // Mu0xce9c
+	0x4E, // Nu0xce9d
 	0xce9e, // Xi
-	0xce9f, // Omicron
-	0xcea0, // Pi
-	0xcea1, // Rho
-	0xcea3, // Sigma
-	0xcea4, // Tau
-	0xcea5, // Upsilon
-	0xcea6, // Phi
-	0xcea7, // Chi
-	0xcea8, // Psi
-	0xcea9, // Omega
+	0x4F, // Omicron0xce9f
+	0xE3, // Pi0xcea0
+	0x50, // Rho0xcea1
+	0xE4, // Sigma 0xcea3
+	0x54, // Tau0xcea4
+	0x59, // Upsilon0xcea5
+	0xE8, // Phi0xcea6
+	0x58, // Chi0xcea7
+	0xcea8, // Psi0xcea8
+	0xEA, // Omega0xcea9
 
-	0xceb1, // alpha
-	0xceb2, // beta
-	0xceb3, // gamma
-	0xceb4, // delta
-	0xceb5, // epsilon
-	0xceb6, // zeta
-	0xceb7, // eta
-	0xceb8, // theta
-	0xceb9, // iota
-	0xceba, // kappa
+	0xE0, // alpha0xceb1
+	0xE1, // beta0xceb2
+	0x72, // gamma0xceb3
+	0xEB, // delta0xceb4
+	0xEE, // epsilon0xceb5
+	0x7A, // zeta0xceb6
+	0x6E, // eta0xceb7
+	0xE9, // theta0xceb8
+	0x69, // iota0xceb9
+	0x6B, // kappa0xceba
 	0xcebb, // lambda
-	0xcebc, // mu
-	0xcebd, // nu
+	0xE6, // mu0xcebc
+	0x76, // nu0xcebd
 	0xcebe, // xi
-	0xcebf, // omicron
-	0xcf80, // pi
-	0xcf81, // rho
-	0xcf83, // sigma
-	0xcf84, // tau
-	0xcf85, // upsilon
-	0xcf86, // phi
-	0xcf87, // chi
+	0x6F, // omicron0xcebf
+	0xE3, // pi0xcf80
+	0x70, // rho0xcf81
+	0xE5, // sigma0xcf83
+	0xE7, // tau0xcf84
+	0x75, // upsilon0xcf85
+	0xED, // phi0xcf86
+	0x78, // chi0xcf87
 	0xcf88, // psi
-	0xcf89, // omega
+	0x77, // omega0xcf89
 
 	0xc4a7, // hbar
 };
@@ -16164,10 +15890,11 @@ fmt_symbol_fragment(char *s, int k)
 	char *t;
 
 	for (i = 0; i < NUM_SYMBOL_NAMES; i++) {
-		t = symbol_name_tab[i];
-		n = (int) strlen(t);
-		if (strncmp(s + k, t, n) == 0)
+		t = (char *)symbol_name_tab[i];
+		n = (int) strlen((const char *)t);
+		if (strncmp(s + k, t, n) == 0){
 			break;
+		}
 	}
 
 	if (i == NUM_SYMBOL_NAMES) {
@@ -16175,7 +15902,7 @@ fmt_symbol_fragment(char *s, int k)
 		return k + 1;
 	}
 
-	c = symbol_unicode_tab[i];
+	c = (int)symbol_unicode_tab[i];
 
 	fmt_roman_char(c);
 
@@ -16229,19 +15956,21 @@ fmt_table(int x, int y, struct atom *p)
 void
 fmt_tensor(struct atom *p)
 {
-	if (p->u.tensor->ndim % 2 == 1)
+	if (p->u.tensor->ndim % 2 == 1){
 		fmt_vector(p); // odd rank
-	else
+	}else{
 		fmt_matrix(p, 0, 0); // even rank
+	}
 }
 
 void
 fmt_term(struct atom *p)
 {
-	if (car(p) == symbol(MULTIPLY))
+	if (car(p) == symbol(MULTIPLY)){
 		fmt_term_nib(p);
-	else
+	}else{
 		fmt_factor(p);
+	}
 }
 
 void
@@ -16256,8 +15985,9 @@ fmt_term_nib(struct atom *p)
 
 	p = cdr(p);
 
-	if (isminusone(car(p)))
+	if (isminusone(car(p))){
 		p = cdr(p); // sign already emitted
+	}
 
 	fmt_factor(car(p));
 
@@ -16302,8 +16032,9 @@ fmt_update_list(int t)
 	int d, h, i, w;
 	struct atom *p1;
 
-	if (tos - t == 1)
+	if (tos - t == 1){
 		return;
+	}
 
 	h = 0;
 	d = 0;
@@ -16522,13 +16253,15 @@ fmt_vector(struct atom *p)
 
 	n = p->u.tensor->ndim;
 
-	for (i = 1; i < n; i++)
+	for (i = 1; i < n; i++){
 		span *= p->u.tensor->dim[i];
+	}
 
 	n = p->u.tensor->dim[0]; // number of rows
 
-	for (i = 0; i < n; i++)
+	for (i = 0; i < n; i++){
 		fmt_matrix(p, 1, i * span);
+	}
 
 	fmt_update_table(n, 1); // n rows, 1 column
 }
@@ -16582,8 +16315,9 @@ fmt_draw(int x, int y, struct atom *p)
 
 		fmt_draw_char(x, y, BDLR);
 
-		for (i = 1; i < w - 1; i++)
+		for (i = 1; i < w - 1; i++){
 			fmt_draw_char(x + i, y, BDLH);
+		}
 
 		fmt_draw_char(x + w - 1, y, BDLL);
 
@@ -16612,8 +16346,9 @@ fmt_draw(int x, int y, struct atom *p)
 void
 fmt_draw_char(int x, int y, int c)
 {
-	if (x >= 0 && x < fmt_ncol && y >= 0 && y < fmt_nrow)
+	if (x >= 0 && x < fmt_ncol && y >= 0 && y < fmt_nrow){
 		fmt_buf[y * fmt_ncol + x] = c;
+	}
 }
 
 void
@@ -16635,8 +16370,9 @@ fmt_draw_ldelim(int x, int y, int h, int d)
 
 	fmt_draw_char(x, y - h + 1, BDLDAR);
 
-	for (i = 1; i < h + d - 1; i++)
+	for (i = 1; i < h + d - 1; i++){
 		fmt_draw_char(x, y - h + 1 + i, BDLV);
+	}
 
 	fmt_draw_char(x, y + d, BDLUAR);
 }
@@ -16648,8 +16384,9 @@ fmt_draw_rdelim(int x, int y, int h, int d)
 
 	fmt_draw_char(x, y - h + 1, BDLDAL);
 
-	for (i = 1; i < h + d - 1; i++)
+	for (i = 1; i < h + d - 1; i++){
 		fmt_draw_char(x, y - h + 1 + i, BDLV);
+	}
 
 	fmt_draw_char(x, y + d, BDLUAL);
 }
@@ -16703,16 +16440,19 @@ void
 fmt_putw(uint32_t w)
 {
 	uint8_t buf[4];
-	if (w == 0)
+	if (w == 0){
 		w = ' ';
+	}
 	buf[0] = w >> 24;
 	buf[1] = w >> 16;
 	buf[2] = w >> 8;
 	buf[3] = w;
-	if (buf[1])
+	if (buf[1]){
 		outbuf_putc(buf[1]);
-	if (buf[2])
+	}
+	if (buf[2]){
 		outbuf_putc(buf[2]);
+	}
 	outbuf_putc(buf[3]);
 }
 // automatic variables not visible to the garbage collector are reclaimed
@@ -16720,7 +16460,7 @@ fmt_putw(uint32_t w)
 void
 gc(void)
 {
-	int i, j;
+	int i;
 	struct atom *p;
 
 	gc_count++;
@@ -16728,9 +16468,9 @@ gc(void)
 
 	// tag everything
 
-	for (i = 0; i < block_count; i++)
-		for (j = 0; j < BLOCKSIZE; j++)
-			mem[i][j].tag = 1;
+	for (i = 0; i < MAXATOMS;i++){
+		mem[i].tag = 1;
+	}
 
 	// untag what's used
 
@@ -16739,8 +16479,9 @@ gc(void)
 	untag(minusone);
 	untag(imaginaryunit);
 
-	for (i = 0; i < tos; i++)
+	for (i = 0; i < tos; i++){
 		untag(stack[i]);
+	}
 
 	for (i = 0; i < 27 * BUCKETSIZE; i++) {
 		untag(symtab[i]);
@@ -16753,23 +16494,20 @@ gc(void)
 	free_list = NULL;
 	free_count = 0;
 
-	for (i = 0; i < block_count; i++)
-		for (j = 0; j < BLOCKSIZE; j++) {
+	for (i = 0; i < MAXATOMS;i++){
 
-			p = mem[i] + j;
-
-			if (p->tag == 0)
-				continue;
-
-			// still tagged so it's unused, put on free list
-
-			switch (p->atomtype) {
+		p = &mem[i];
+		if (p->tag == 0){
+			continue;
+		}
+		// still tagged so it's unused, put on free list
+		switch (p->atomtype) {
 			case KSYM:
-				m_free(p->u.ksym.name);
+				e_free(p->u.ksym.name);
 				ksym_count--;
 				break;
 			case USYM:
-				m_free(p->u.usym.name);
+				e_free(p->u.usym.name);
 				usym_count--;
 				break;
 			case RATIONAL:
@@ -16778,23 +16516,21 @@ gc(void)
 				break;
 			case STR:
 				if (p->u.str)
-					m_free(p->u.str);
+					e_free(p->u.str);
 				string_count--;
 				break;
 			case TENSOR:
-				m_free(p->u.tensor);
+				e_free(p->u.tensor);
 				tensor_count--;
 				break;
 			default:
 				break; // FREEATOM, CONS, or DOUBLE
-			}
-
-			p->atomtype = FREEATOM;
-			p->u.next = free_list;
-
-			free_list = p;
-			free_count++;
 		}
+		p->atomtype = FREEATOM;
+		p->u.next = free_list;
+		free_list = p;
+		free_count++;
+	}
 }
 
 void
@@ -16802,36 +16538,41 @@ untag(struct atom *p)
 {
 	int i;
 
-	if (p == NULL)
+	if (p == NULL){
 		return;
+	}
 
 	while (iscons(p)) {
-		if (p->tag == 0)
+		if (p->tag == 0){
 			return;
+		}
 		p->tag = 0;
 		untag(p->u.cons.car);
 		p = p->u.cons.cdr;
 	}
 
-	if (p->tag == 0)
+	if (p->tag == 0){
 		return;
+	}
 
 	p->tag = 0;
 
-	if (istensor(p))
-		for (i = 0; i < p->u.tensor->nelem; i++)
+	if (istensor(p)){
+		for (i = 0; i < p->u.tensor->nelem; i++){
 			untag(p->u.tensor->elem[i]);
+		}
+	}
 }
-struct atom *mem[MAXBLOCKS]; // an array of pointers
+struct atom *mem; // an array of pointers
 struct atom *free_list;
 
 int tos; // top of stack
 
-struct atom *stack[STACKSIZE];
+struct atom **stack;
 
-struct atom *symtab[27 * BUCKETSIZE];
-struct atom *binding[27 * BUCKETSIZE];
-struct atom *usrfunc[27 * BUCKETSIZE];
+struct atom **symtab;
+struct atom **binding;
+struct atom **usrfunc;
 
 struct atom *zero;
 struct atom *one;
@@ -16842,11 +16583,10 @@ int eval_level;
 int gc_level;
 int expanding;
 int drawing;
+int nonstop;
 int interrupt;
-int shuntflag;
-int errorflag;
-int breakflag;
-jmp_buf jmpbuf;
+jmp_buf jmpbuf0;
+jmp_buf jmpbuf1;
 char *trace1;
 char *trace2;
 
@@ -16867,51 +16607,44 @@ char strbuf[STRBUFLEN];
 char *outbuf;
 int outbuf_index;
 int outbuf_length;
-// int
-// main(int argc, char *argv[])
-// {
-// 	int i;
-// 	for (i = 1; i < argc; i++)
-// 		run_infile(argv[i]);
-// 	if (isatty(fileno(stdout)))
-// 		run_stdin();
-// }
-
+/*
+int
+main(int argc, char *argv[])
+{
+	int i;
+	for (i = 1; i < argc; i++)
+		run_infile(argv[i]);
+	if (isatty(fileno(stdout)))
+		run_stdin();
+}
+*/
 void
 run_infile(char *infile)
 {
 	char *buf;
 	buf = read_file(infile);
 	if (buf == NULL) {
-		fprintf(stderr, "cannot read %s\n", infile);
-		exit(1);
+		//fprintf(stderr, "cannot read %s\n", infile);
+		mp_printf(&mp_plat_print, "cannot read %s\n", infile);
+		//exit(1);
+		stopf("cannot read file");
 	}
 	run(buf);
-	m_free(buf);
+	e_free(buf);
 }
 
 void
 run_stdin(void)
 {
-	mp_printf(&mp_plat_print, "Eigenmath start %d\n", __LINE__);
-	// static char inbuf[1000];
-	vstr_t* vstr_inbuf = vstr_new(1);
-
-	// for (;;) {
-		// fputs("? ", stdout);
-		// fflush(stdout);
-		// fgets(inbuf, sizeof inbuf, stdin);
-		int res = readline(vstr_inbuf,"eigenmath> ");
-		mp_printf(&mp_plat_print, "Eigenmath run:\n");
-		mp_printf(&mp_plat_print, "res=%d\n", res);
-		mp_printf(&mp_plat_print, "%s\n", vstr_inbuf->buf);
-		run(vstr_inbuf->buf);
-	// }
+	static char inbuf[1000];
+	for (;;) {
+		fputs("? ", stdout);
+		fflush(stdout);
+		fgets(inbuf, sizeof inbuf, stdin);
+		run(inbuf);
+	}
 }
 
-
-
-
 void
 display(void)
 {
@@ -16921,8 +16654,21 @@ display(void)
 void
 printbuf(char *s, int color)
 {
-	// fputs(s, stdout);
-	mp_printf(&mp_plat_print, "%s", s);
+	//fputs(s, stdout);
+
+	switch (color) {
+		case 0:
+			mp_printf(&mp_plat_print, "\x1b[37;40m%s\x1b[0m", s);//black
+			break;
+		case 1:
+			mp_printf(&mp_plat_print, "\x1b[37;44m%s\x1b[0m", s); // blue
+			break;
+		case 2:
+			mp_printf(&mp_plat_print, "\x1b[37;41m%s\x1b[0m", s); // red
+			break;
+	}
+
+
 }
 
 void
@@ -16936,7 +16682,8 @@ void
 eval_exit(struct atom *p1)
 {
 	(void) p1; // silence compiler
-	exit(0);
+	//exit(0);
+	longjmp(jmpbuf0, 1);
 }
 void
 numden(void)
@@ -16970,8 +16717,9 @@ numden_find_divisor(struct atom *p)
 	if (car(p) == symbol(ADD)) {
 		p = cdr(p);
 		while (iscons(p)) {
-			if (numden_find_divisor_term(car(p)))
+			if (numden_find_divisor_term(car(p))){
 				return 1;
+			}
 			p = cdr(p);
 		}
 		return 0;
@@ -16986,8 +16734,9 @@ numden_find_divisor_term(struct atom *p)
 	if (car(p) == symbol(MULTIPLY)) {
 		p = cdr(p);
 		while (iscons(p)) {
-			if (numden_find_divisor_factor(car(p)))
+			if (numden_find_divisor_factor(car(p))){
 				return 1;
+			}
 			p = cdr(p);
 		}
 		return 0;
@@ -17005,9 +16754,9 @@ numden_find_divisor_factor(struct atom *p)
 	}
 
 	if (car(p) == symbol(POWER) && isnegativeterm(caddr(p))) {
-		if (isminusone(caddr(p)))
+		if (isminusone(caddr(p))){
 			push(cadr(p));
-		else {
+		}else {
 			push_symbol(POWER);
 			push(cadr(p));
 			push(caddr(p));
@@ -17055,30 +16804,34 @@ outbuf_init(void)
 	outbuf_puts(""); // init outbuf as empty string
 }
 
-void
-outbuf_puts(char *s)
+
+void outbuf_puts(char *s)
 {
-	int len, m;
+    int len, m;
 
-	len = (int) strlen(s);
+    len = (int) strlen((const char *)s);
 
-	// Let outbuf_index + len == 1000
+    // Make sure there is enough room for new string + '\0' terminator
+    m = 1000 * ((outbuf_index + len + 1) / 1000 + 1); // +1 for '\0'
 
-	// Then m == 2000 hence there is always room for the terminator '\0'
+    if (m > outbuf_length) {
+        outbuf = e_realloc(outbuf, m);
+        if (outbuf == NULL) {
+            stopf("outbuf_puts: realloc failed");
+        }
+        outbuf_length = m;
+    }
 
-	m = 1000 * ((outbuf_index + len) / 1000 + 1); // m is a multiple of 1000
+    // Copy the string
+    memcpy(outbuf + outbuf_index, s, len);
+    outbuf_index += len;
 
-	if (m > outbuf_length) {
-		outbuf = realloc(outbuf, m);
-		if (outbuf == NULL)
-			exit(1);
-		outbuf_length = m;
-	}
-
-	strcpy(outbuf + outbuf_index, s);
-	outbuf_index += len;
+    // Always null-terminate
+    outbuf[outbuf_index] = '\0';
 }
 
+
+
 void
 outbuf_putc(int c)
 {
@@ -17091,9 +16844,11 @@ outbuf_putc(int c)
 	m = 1000 * ((outbuf_index + 1) / 1000 + 1); // m is a multiple of 1000
 
 	if (m > outbuf_length) {
-		outbuf = realloc(outbuf, m);
-		if (outbuf == NULL)
-			exit(1);
+		outbuf = e_realloc(outbuf, m);
+		if (outbuf == NULL){
+			//exit(1);
+			stopf("outbuf_putc: realloc failed");
+		}
 		outbuf_length = m;
 	}
 
@@ -17104,14 +16859,18 @@ int
 iszero(struct atom *p)
 {
 	int i;
-	if (isrational(p))
+	if (isrational(p)){
 		return MZERO(p->u.q.a);
-	if (isdouble(p))
+	}
+	if (isdouble(p)){
 		return p->u.d == 0.0;
+	}
 	if (istensor(p)) {
-		for (i = 0; i < p->u.tensor->nelem; i++)
-			if (!iszero(p->u.tensor->elem[i]))
+		for (i = 0; i < p->u.tensor->nelem; i++){
+			if (!iszero(p->u.tensor->elem[i])){
 				return 0;
+			}
+		}
 		return 1;
 	}
 	return 0;
@@ -17131,12 +16890,14 @@ isequalq(struct atom *p, int a, int b)
 		if (a < 0) {
 			sign = MMINUS;
 			a = -a;
-		} else
+		} else{
 			sign = MPLUS;
+		}
 		return p->sign == sign && MEQUAL(p->u.q.a, a) && MEQUAL(p->u.q.b, b);
 	}
-	if (isdouble(p))
+	if (isdouble(p)){
 		return p->u.d == (double) a / b;
+	}
 	return 0;
 }
 
@@ -17191,24 +16952,27 @@ isnegativeterm(struct atom *p)
 int
 isnegativenumber(struct atom *p)
 {
-	if (isrational(p))
+	if (isrational(p)){
 		return p->sign == MMINUS;
-	else if (isdouble(p))
+	}else if (isdouble(p)){
 		return p->u.d < 0.0;
-	else
+	}else{
 		return 0;
+	}
 }
 
 int
 isimaginaryterm(struct atom *p)
 {
-	if (isimaginaryfactor(p))
+	if (isimaginaryfactor(p)){
 		return 1;
+	}
 	if (car(p) == symbol(MULTIPLY)) {
 		p = cdr(p);
 		while (iscons(p)) {
-			if (isimaginaryfactor(car(p)))
+			if (isimaginaryfactor(car(p))){
 				return 1;
+			}
 			p = cdr(p);
 		}
 	}
@@ -17258,34 +17022,41 @@ isdoublez(struct atom *p)
 {
 	if (car(p) == symbol(ADD)) {
 
-		if (lengthf(p) != 3)
+		if (lengthf(p) != 3){
 			return 0;
+		}
 
-		if (!isdouble(cadr(p))) // x
+		if (!isdouble(cadr(p))) {// x
 			return 0;
+		}
 
 		p = caddr(p);
 	}
 
-	if (car(p) != symbol(MULTIPLY))
+	if (car(p) != symbol(MULTIPLY)){
 		return 0;
+	}
 
-	if (lengthf(p) != 3)
+	if (lengthf(p) != 3){
 		return 0;
+	}
 
-	if (!isdouble(cadr(p))) // y
+	if (!isdouble(cadr(p))) {// y
 		return 0;
-
+	}
 	p = caddr(p);
 
-	if (car(p) != symbol(POWER))
+	if (car(p) != symbol(POWER)){
 		return 0;
+	}	
 
-	if (!isminusone(cadr(p)))
+	if (!isminusone(cadr(p))){
 		return 0;
+	}
 
-	if (!isequalq(caddr(p), 1, 2))
+	if (!isequalq(caddr(p), 1, 2)){
 		return 0;
+	}
 
 	return 1;
 }
@@ -17293,23 +17064,25 @@ isdoublez(struct atom *p)
 int
 isdenominator(struct atom *p)
 {
-	if (car(p) == symbol(POWER) && isnegativenumber(caddr(p)))
+	if (car(p) == symbol(POWER) && isnegativenumber(caddr(p))){
 		return 1;
-	else if (isrational(p) && !MEQUAL(p->u.q.b, 1))
+	}else if (isrational(p) && !MEQUAL(p->u.q.b, 1)){
 		return 1;
-	else
+	}else{
 		return 0;
+	}
 }
 
 int
 isnumerator(struct atom *p)
 {
-	if (car(p) == symbol(POWER) && isnegativenumber(caddr(p)))
+	if (car(p) == symbol(POWER) && isnegativenumber(caddr(p))){
 		return 0;
-	else if (isrational(p) && MEQUAL(p->u.q.a, 1))
+	} else if (isrational(p) && MEQUAL(p->u.q.a, 1)){
 		return 0;
-	else
+	}else{
 		return 1;
+	}
 }
 
 int
@@ -17318,8 +17091,9 @@ hasdouble(struct atom *p)
 	if (iscons(p)) {
 		p = cdr(p);
 		while (iscons(p)) {
-			if (hasdouble(car(p)))
+			if (hasdouble(car(p))){
 				return 1;
+			}
 			p = cdr(p);
 		}
 		return 0;
@@ -17333,8 +17107,9 @@ isdenormalpolar(struct atom *p)
 	if (car(p) == symbol(ADD)) {
 		p = cdr(p);
 		while (iscons(p)) {
-			if (isdenormalpolarterm(car(p)))
+			if (isdenormalpolarterm(car(p))){
 				return 1;
+			}
 			p = cdr(p);
 		}
 		return 0;
@@ -17348,28 +17123,31 @@ isdenormalpolar(struct atom *p)
 int
 isdenormalpolarterm(struct atom *p)
 {
-	if (car(p) != symbol(MULTIPLY))
+	if (car(p) != symbol(MULTIPLY)){
 		return 0;
+	}
 
-	if (lengthf(p) == 3 && isimaginaryunit(cadr(p)) && caddr(p) == symbol(PI))
+	if (lengthf(p) == 3 && isimaginaryunit(cadr(p)) && caddr(p) == symbol(PI)){
 		return 1; // exp(i pi)
+		}
 
-	if (lengthf(p) != 4 || !isnum(cadr(p)) || !isimaginaryunit(caddr(p)) || cadddr(p) != symbol(PI))
+	if (lengthf(p) != 4 || !isnum(cadr(p)) || !isimaginaryunit(caddr(p)) || cadddr(p) != symbol(PI)){
 		return 0;
+	}
 
 	p = cadr(p); // p = coeff of term
 
-	if (isnegativenumber(p))
+	if (isnegativenumber(p)){
 		return 1; // p < 0
-
+	}
 	push(p);
 	push_rational(-1, 2);
 	add();
 	p = pop();
 
-	if (!isnegativenumber(p))
+	if (!isnegativenumber(p)){
 		return 1; // p >= 1/2
-
+	}
 	return 0;
 }
 
@@ -17382,55 +17160,31 @@ issquarematrix(struct atom *p)
 int
 issmallinteger(struct atom *p)
 {
-	if (isinteger(p))
+	if (isinteger(p)){
 		return MLENGTH(p->u.q.a) == 1 && p->u.q.a[0] <= 0x7fffffff;
+	}
 
-	if (isdouble(p))
+	if (isdouble(p)){
 		return p->u.d == floor(p->u.d) && fabs(p->u.d) <= 0x7fffffff;
-
-	return 0;
-}
-
-// does f depend on x?
-
-int
-dependent(struct atom *f, struct atom *x)
-{
-	if (equal(f, x))
-		return 1;
-
-	// a user function with no arguments?
-
-	if (isusersymbol(car(f)) && lengthf(f) == 1)
-		return 1;
-
-	while (iscons(f)) {
-		if (dependent(car(f), x))
-			return 1;
-		f = cdr(f);
 	}
 
 	return 0;
 }
-
 void
 run(char *buf)
 {
-	if (setjmp(jmpbuf))
+	if (setjmp(jmpbuf0)){
 		return;
-
+	}
 	tos = 0;
 	interrupt = 0;
 	eval_level = 0;
 	gc_level = 0;
 	expanding = 1;
 	drawing = 0;
-	shuntflag = 0;
-	errorflag = 0;
-	breakflag = 0;
-	mp_printf(&mp_plat_print, "run() has param: %s\n", buf);
+	nonstop = 0;
+
 	if (zero == NULL) {
-		srand((unsigned) time(NULL));
 		init_symbol_table();
 		push_bignum(MPLUS, mint(0), mint(1));
 		zero = pop();
@@ -17443,10 +17197,10 @@ run(char *buf)
 		push_rational(1, 2);
 		list(3);
 		imaginaryunit = pop();
-		// run_init_script();
+		run_init_script();
 	}
-	mp_printf(&mp_plat_print, "start to run run_buf()\n");
-	// run_buf(buf);
+
+	run_buf(buf);
 }
 
 void
@@ -17475,8 +17229,9 @@ run_buf(char *buf)
 		dupl();
 		p1 = pop();
 
-		if (p1 != symbol(NIL))
+		if (p1 != symbol(NIL)){
 			set_symbol(symbol(LAST), p1, symbol(NIL));
+		}
 
 		print_result();
 	}
@@ -17493,8 +17248,9 @@ scan_input(char *s)
 	s = scan(s);
 	trace2 = s;
 	p1 = get_binding(symbol(TRACE));
-	if (p1 != symbol(TRACE) && !iszero(p1))
+	if (p1 != symbol(TRACE) && !iszero(p1)){
 		print_trace(BLUE);
+	}
 	return s;
 }
 
@@ -17511,12 +17267,17 @@ print_trace(int color)
 		c = *s++;
 		outbuf_putc(c);
 	}
-	if (c != '\n')
+	if (c != '\n'){
 		outbuf_putc('\n');
-	printbuf(outbuf, color);
+	}
+	if (noprint == false){
+		printbuf(outbuf, color);
+	}
+	
 }
 
-char *init_script =
+const char *init_script =
+"tty = 1\n"
 "i = sqrt(-1)\n"
 "grad(f) = d(f,(x,y,z))\n"
 "cross(a,b) = (dot(a[2],b[3])-dot(a[3],b[2]),dot(a[3],b[1])-dot(a[1],b[3]),dot(a[1],b[2])-dot(a[2],b[1]))\n"
@@ -17532,16 +17293,32 @@ char *init_script =
 void
 run_init_script(void)
 {
-	run_buf(init_script);
+	run_buf((char *)init_script);
 }
 
 void
 stopf(char *s)
 {
+	if (nonstop){
+		longjmp(jmpbuf1, 1);
+	}
 	print_trace(RED);
 	snprintf(strbuf, STRBUFLEN, "Stop: %s\n", s);
-	printbuf(strbuf, RED);
-	longjmp(jmpbuf, 1);
+	if (noprint == false){
+		printbuf(strbuf, RED);
+	}else{
+		//todo:push into outbuf
+		outbuf_puts(strbuf);
+	}
+	
+	longjmp(jmpbuf0, 1);
+}
+
+void
+exitf(char *s)
+{
+	nonstop = 0;
+	stopf(s);
 }
 // token_str and scan_str are pointers to the input string, for example
 //
@@ -17606,11 +17383,13 @@ scan_nib(char *s)
 	scan_str = s;
 	scan_level = 0;
 	get_token_skip_newlines();
-	if (token == T_END)
+	if (token == T_END){
 		return NULL;
+	}
 	scan_stmt();
-	if (token != T_NEWLINE && token != T_END)
+	if (token != T_NEWLINE && token != T_END){
 		scan_error("syntax err"); // unexpected token, for example, 1:2
+	}
 	return scan_str;
 }
 
@@ -17632,23 +17411,23 @@ scan_comparison(void)
 {
 	scan_expression();
 	switch (token) {
-	case T_EQ:
-		push_symbol(TESTEQ); // ==
-		break;
-	case T_LTEQ:
-		push_symbol(TESTLE);
-		break;
-	case T_GTEQ:
-		push_symbol(TESTGE);
-		break;
-	case '<':
-		push_symbol(TESTLT);
-		break;
-	case '>':
-		push_symbol(TESTGT);
-		break;
-	default:
-		return;
+		case T_EQ:
+			push_symbol(TESTEQ); // ==
+			break;
+		case T_LTEQ:
+			push_symbol(TESTLE);
+			break;
+		case T_GTEQ:
+			push_symbol(TESTGE);
+			break;
+		case '<':
+			push_symbol(TESTLT);
+			break;
+		case '>':
+			push_symbol(TESTGT);
+			break;
+		default:
+			return;
 	}
 	swap();
 	get_token_skip_newlines(); // get token after rel op
@@ -17661,17 +17440,20 @@ scan_expression(void)
 {
 	int h = tos, t;
 	t = token;
-	if (token == '+' || token == '-')
+	if (token == '+' || token == '-'){
 		get_token_skip_newlines();
+	}
 	scan_term();
-	if (t == '-')
+	if (t == '-'){
 		static_negate();
+	}
 	while (token == '+' || token == '-') {
 		t = token;
 		get_token_skip_newlines(); // get token after '+' or '-'
 		scan_term();
-		if (t == '-')
+		if (t == '-'){
 			static_negate();
+		}
 	}
 	if (tos - h > 1) {
 		list(tos - h);
@@ -17685,17 +17467,17 @@ int
 another_factor_pending(void)
 {
 	switch (token) {
-	case '*':
-	case '/':
-	case '(':
-	case T_SYMBOL:
-	case T_FUNCTION:
-	case T_INTEGER:
-	case T_DOUBLE:
-	case T_STRING:
-		return 1;
-	default:
-		break;
+		case '*':
+		case '/':
+		case '(':
+		case T_SYMBOL:
+		case T_FUNCTION:
+		case T_INTEGER:
+		case T_DOUBLE:
+		case T_STRING:
+			return 1;
+		default:
+			break;
 	}
 	return 0;
 }
@@ -17711,13 +17493,15 @@ scan_term(void)
 
 		t = token;
 
-		if (token == '*' || token == '/')
+		if (token == '*' || token == '/'){
 			get_token_skip_newlines();
+		}
 
 		scan_power();
 
-		if (t == '/')
+		if (t == '/'){
 			static_reciprocate();
+		}
 	}
 
 	if (tos - h > 1) {
@@ -17791,8 +17575,9 @@ scan_factor(void)
 			get_token(); // get token after ','
 			scan_expression();
 		}
-		if (token != ']')
+		if (token != ']'){
 			scan_error("expected ']'");
+		}
 		scan_level--;
 		get_token(); // get token after ']'
 		list(tos - h);
@@ -17809,7 +17594,7 @@ scan_factor(void)
 void
 scan_symbol(void)
 {
-	if (scan_mode && strlen(token_buf) == 1)
+	if (scan_mode && strlen((const char *)token_buf) == 1){
 		switch (token_buf[0]) {
 		case 'a':
 			push_symbol(SA);
@@ -17824,8 +17609,9 @@ scan_symbol(void)
 			push(lookup(token_buf));
 			break;
 		}
-	else
+	}else{
 		push(lookup(token_buf));
+	}
 	get_token();
 }
 
@@ -17855,8 +17641,9 @@ scan_function_call(void)
 		get_token(); // get token after ','
 		scan_stmt();
 	}
-	if (token != ')')
+	if (token != ')'){
 		scan_error("expected ')'");
+	}
 	scan_level--;
 	get_token(); // get token after ')'
 	list(tos - h);
@@ -17897,16 +17684,19 @@ scan_subexpr(void)
 		get_token(); // get token after ','
 		scan_stmt();
 	}
-	if (token != ')')
+	if (token != ')'){
 		scan_error("expected ')'");
+	}
 	scan_level--;
 	get_token(); // get token after ')'
 	n = tos - h;
-	if (n < 2)
+	if (n < 2){
 		return;
+	}
 	p = alloc_vector(n);
-	for (i = 0; i < n; i++)
+	for (i = 0; i < n; i++){
 		p->u.tensor->elem[i] = stack[h + i];
+	}
 	tos = h;
 	push(p);
 }
@@ -17923,9 +17713,11 @@ void
 get_token(void)
 {
 	get_token_nib();
-	if (scan_level)
-		while (token == T_NEWLINE)
+	if (scan_level){
+		while (token == T_NEWLINE){
 			get_token_nib();
+		}
+	}
 }
 
 void
@@ -17933,9 +17725,9 @@ get_token_nib(void)
 {
 	// skip spaces
 
-	while (*scan_str != '\0' && *scan_str != '\n' && *scan_str != '\r' && (*scan_str < 33 || *scan_str > 126))
+	while (*scan_str != '\0' && *scan_str != '\n' && *scan_str != '\r' && (*scan_str < 33 || *scan_str > 126)){
 		scan_str++;
-
+	}
 	token_str = scan_str;
 
 	// end of input?
@@ -17956,41 +17748,48 @@ get_token_nib(void)
 	// comment?
 
 	if (*scan_str == '#' || (scan_str[0] == '-' && scan_str[1] == '-')) {
-		while (*scan_str && *scan_str != '\n' && *scan_str != '\r')
+		while (*scan_str && *scan_str != '\n' && *scan_str != '\r'){
 			scan_str++;
-		if (*scan_str)
+		}
+		if (*scan_str){
 			scan_str++;
+		}
 		token = T_NEWLINE;
 		return;
 	}
 
 	// number?
 
-	if (isdigit((unsigned char)*scan_str) || *scan_str == '.') {
-		while (isdigit((unsigned char)*scan_str))
+	if (isdigit((int)*scan_str) || *scan_str == '.') {
+		while (isdigit((int)*scan_str)){
 			scan_str++;
+		}
 		if (*scan_str == '.') {
 			scan_str++;
-			while (isdigit((unsigned char)*scan_str))
+			while (isdigit((int)*scan_str)){
 				scan_str++;
-			if (token_str + 1 == scan_str)
+			}
+			if (token_str + 1 == scan_str){
 				scan_error("expected decimal digit"); // only a decimal point
+			}
 			token = T_DOUBLE;
-		} else
+		} else{
 			token = T_INTEGER;
+		}
 		update_token_buf(token_str, scan_str);
 		return;
 	}
 
 	// symbol?
 
-	if (isalpha((unsigned char)*scan_str)) {
-		while (isalnum((unsigned char)*scan_str))
+	if (isalpha((int)*scan_str)) {
+		while (isalnum((int)*scan_str))
 			scan_str++;
-		if (*scan_str == '(')
+		if (*scan_str == '('){
 			token = T_FUNCTION;
-		else
+		}else{
 			token = T_SYMBOL;
+		}
 		update_token_buf(token_str, scan_str);
 		return;
 	}
@@ -18001,8 +17800,9 @@ get_token_nib(void)
 		scan_str++;
 		while (*scan_str && *scan_str != '"' && *scan_str != '\n' && *scan_str != '\r')
 			scan_str++;
-		if (*scan_str != '"')
+		if (*scan_str != '"'){
 			scan_error("runaway string");
+		}	
 		scan_str++;
 		token = T_STRING;
 		update_token_buf(token_str + 1, scan_str - 1); // don't include quote chars
@@ -18048,8 +17848,9 @@ update_token_buf(char *a, char *b)
 	m = 1000 * (n / 1000 + 1); // m is a multiple of 1000
 
 	if (m > token_buf_len) {
-		if (token_buf)
-			m_free(token_buf);
+		if (token_buf){
+			e_free(token_buf);
+		}
 		token_buf = alloc_mem(m);
 		token_buf_len = m;
 	}
@@ -18125,8 +17926,9 @@ static_reciprocate(void)
 		return;
 	}
 
-	if (!isrational(p1) || !isequaln(p1, 1))
+	if (!isrational(p1) || !isequaln(p1, 1)){
 		push(p1); // p1 != 1
+	}
 
 	if (isnum(p2)) {
 		push(p2);
@@ -18151,18 +17953,21 @@ static_reciprocate(void)
 void
 push(struct atom *p)
 {
-	if (tos < 0 || tos >= STACKSIZE)
-		stopf("stack error");
+	if (tos < 0 || tos >= STACKSIZE){
+		exitf("stack error, circular definition?");
+	}
 	stack[tos++] = p;
-	if (tos > max_tos)
+	if (tos > max_tos){
 		max_tos = tos;
+	}
 }
 
 struct atom *
 pop(void)
 {
-	if (tos < 1 || tos > STACKSIZE)
-		stopf("stack error");
+	if (tos < 1 || tos > STACKSIZE){
+		exitf("stack error");
+	}
 	return stack[--tos];
 }
 
@@ -18228,10 +18033,11 @@ push_integer(int n)
 void
 push_rational(int a, int b)
 {
-	if (a < 0)
+	if (a < 0){
 		push_bignum(MMINUS, mint(-a), mint(b));
-	else
+	}else{
 		push_bignum(MPLUS, mint(a), mint(b));
+	}
 }
 
 void
@@ -18273,9 +18079,9 @@ pop_integer(void)
 		n = p->u.q.a[0];
 		if (isnegativenumber(p))
 			n = -n;
-	} else
+	} else{
 		n = (int) p->u.d;
-
+	}
 	return n;
 }
 
@@ -18297,12 +18103,12 @@ pop_double(void)
 
 	p = pop();
 
-	if (!isnum(p))
+	if (!isnum(p)){
 		stopf("number expected");
-
-	if (isdouble(p))
+	}
+	if (isdouble(p)){
 		d = p->u.d;
-	else {
+	}else {
 		a = mfloat(p->u.q.a);
 		b = mfloat(p->u.q.b);
 		d = a / b;
@@ -18318,10 +18124,13 @@ push_string(char *s)
 {
 	struct atom *p;
 	p = alloc_str();
-	s = strdup(s);
-	if (s == NULL)
-		exit(1);
-	p->u.str = s;
+	//s = strdup(s);
+	char *ns = e_malloc(strlen((const char *)s) + 1);
+	if (ns == NULL){
+		stopf("push_string memory alloc error");
+	}	
+	memcpy(ns, s, strlen((const char *)s) + 1);
+	p->u.str = ns;
 	push(p);
 }
 
@@ -18332,10 +18141,12 @@ slice(int h, int n)
 {
 	int i, m;
 	m = tos - h - n;
-	if (m < 0)
+	if (m < 0){
 		stopf("stack slice error");
-	for (i = 0; i < m; i++)
+	}	
+	for (i = 0; i < m; i++){
 		stack[h + i] = stack[h + n + i];
+	}	
 	tos -= n;
 }
 // symbol lookup, create symbol if not found
@@ -18346,30 +18157,37 @@ lookup(char *s)
 	int i, k;
 	struct atom *p;
 
-	if (isupper((unsigned char)*s))
+	if (isupper((int)*s)){
 		k = BUCKETSIZE * (*s - 'A');
-	else if (islower((unsigned char)*s))
+	}else if (islower((int)*s)){
 		k = BUCKETSIZE * (*s - 'a');
-	else
+	}else{
 		k = BUCKETSIZE * 26;
-
+	}
 	for (i = 0; i < BUCKETSIZE; i++) {
 		p = symtab[k + i];
-		if (p == NULL)
+		if (p == NULL){
 			break;
-		if (strcmp(s, printname(p)) == 0)
+		}
+		if (strcmp(s, printname(p)) == 0){
 			return p;
+		}	
 	}
 
-	if (i == BUCKETSIZE)
+	if (i == BUCKETSIZE){
 		stopf("symbol table full");
-
+	}
 	p = alloc_atom();
-	s = strdup(s);
-	if (s == NULL)
-		exit(1);
+	//s = strdup(s);
+	char *ns = e_malloc(strlen((const char *)s) + 1);
+	if (ns == NULL){
+		stopf("lookup memory alloc error");
+	}
+	memcpy(ns, s, strlen((const char *)s) + 1);
+	//if (s == NULL)
+	//	exit(1);
 	p->atomtype = USYM;
-	p->u.usym.name = s;
+	p->u.usym.name = ns;
 	p->u.usym.index = k + i;
 	symtab[k + i] = p;
 	usym_count++;
@@ -18380,19 +18198,22 @@ lookup(char *s)
 char *
 printname(struct atom *p)
 {
-	if (iskeyword(p))
+	if (iskeyword(p)){
 		return p->u.ksym.name;
-	else if (isusersymbol(p))
+	}else if (isusersymbol(p)){
 		return p->u.usym.name;
-	else
+	}else{
 		return "?";
+	}		
 }
 
 void
 set_symbol(struct atom *p1, struct atom *p2, struct atom *p3)
 {
-	if (!isusersymbol(p1))
+	if (!isusersymbol(p1)){
 		stopf("symbol error");
+	}
+		
 	binding[p1->u.usym.index] = p2;
 	usrfunc[p1->u.usym.index] = p3;
 }
@@ -18401,30 +18222,35 @@ struct atom *
 get_binding(struct atom *p1)
 {
 	struct atom *p2;
-	if (!isusersymbol(p1))
+	if (!isusersymbol(p1)){
 		stopf("symbol error");
+	}
 	p2 = binding[p1->u.usym.index];
-	if (p2 == NULL || p2 == symbol(NIL))
+	if (p2 == NULL || p2 == symbol(NIL)){
 		p2 = p1; // symbol binds to itself
+	}
 	return p2;
 }
 
 struct atom *
 get_usrfunc(struct atom *p)
 {
-	if (!isusersymbol(p))
+	if (!isusersymbol(p)){
 		stopf("symbol error");
+	}
 	p = usrfunc[p->u.usym.index];
-	if (p == NULL)
+	if (p == NULL){
 		p = symbol(NIL);
+	}
 	return p;
 }
 
-struct se {
-	char *str;
+typedef struct{
+	const char *str;
 	int index;
 	void (*func)(struct atom *);
-} stab[] = {
+} se; 
+const se stab[] = {
 
 	{ "abs",		ABS,		eval_abs		},
 	{ "adj",		ADJ,		eval_adj		},
@@ -18438,7 +18264,6 @@ struct se {
 	{ "arg",		ARG,		eval_arg		},
 
 	{ "binding",		BINDING,	eval_binding		},
-	{ "break",		BREAK,		eval_break		},
 
 	{ "C",			C_UPPER,	NULL			},
 	{ "c",			C_LOWER,	NULL			},
@@ -18501,9 +18326,7 @@ struct se {
 	{ "kronecker",		KRONECKER,	eval_kronecker		},
 
 	{ "last",		LAST,		NULL			},
-	{ "lgamma",		LGAMMA,		eval_lgamma		},
 	{ "log",		LOG,		eval_log		},
-	{ "loop",		LOOP,		eval_loop		},
 
 	{ "mag",		MAG,		eval_mag		},
 	{ "minor",		MINOR,		eval_minor		},
@@ -18534,7 +18357,6 @@ struct se {
 
 	{ "R",			R_UPPER,	NULL			},
 	{ "r",			R_LOWER,	NULL			},
-	{ "rand",		RAND,		eval_rand		},
 	{ "rank",		RANK,		eval_rank		},
 	{ "rationalize",	RATIONALIZE,	eval_rationalize	},
 	{ "real",		REAL,		eval_real		},
@@ -18565,7 +18387,6 @@ struct se {
 	{ "testgt",		TESTGT,		eval_testgt		},
 	{ "testle",		TESTLE,		eval_testle		},
 	{ "testlt",		TESTLT,		eval_testlt		},
-	{ "tgamma",		TGAMMA,		eval_tgamma		},
 	{ "trace",		TRACE,		NULL			},
 	{ "transpose",		TRANSPOSE,	eval_transpose		},
 	{ "tty",		TTY,		NULL			},
@@ -18599,15 +18420,15 @@ struct se {
 	{ "$a",			SA,		NULL			},
 	{ "$b",			SB,		NULL			},
 	{ "$x",			SX,		NULL			},
-	{ "$arg1",		ARG1,		NULL			},
-	{ "$arg2",		ARG2,		NULL			},
-	{ "$arg3",		ARG3,		NULL			},
-	{ "$arg4",		ARG4,		NULL			},
-	{ "$arg5",		ARG5,		NULL			},
-	{ "$arg6",		ARG6,		NULL			},
-	{ "$arg7",		ARG7,		NULL			},
-	{ "$arg8",		ARG8,		NULL			},
-	{ "$arg9",		ARG9,		NULL			},
+	{ "$1",			ARG1,		NULL			},
+	{ "$2",			ARG2,		NULL			},
+	{ "$3",			ARG3,		NULL			},
+	{ "$4",			ARG4,		NULL			},
+	{ "$5",			ARG5,		NULL			},
+	{ "$6",			ARG6,		NULL			},
+	{ "$7",			ARG7,		NULL			},
+	{ "$8",			ARG8,		NULL			},
+	{ "$9",			ARG9,		NULL			},
 };
 
 void
@@ -18623,13 +18444,17 @@ init_symbol_table(void)
 		usrfunc[i] = NULL;
 	}
 
-	n = sizeof stab / sizeof (struct se);
+	n = sizeof stab / sizeof (se);
 
 	for (i = 0; i < n; i++) {
 		p = alloc_atom();
-		s = strdup(stab[i].str);
-		if (s == NULL)
-			exit(1);
+		s = e_malloc(strlen((const char *)(stab[i].str)) + 1);
+		if (s == NULL){
+			stopf("symbol table init error");
+		}
+			
+		memcpy(s, stab[i].str, strlen((const char *)(stab[i].str)) + 1);	
+		
 		if (stab[i].func) {
 			p->atomtype = KSYM;
 			p->u.ksym.name = s;
@@ -18644,182 +18469,3 @@ init_symbol_table(void)
 		symtab[stab[i].index] = p;
 	}
 }
-
-
-// This is the function which will be called from Python as cexample.add_ints(a, b).
-static mp_obj_t example_add_ints(mp_obj_t a_obj, mp_obj_t b_obj) {
-    // Extract the ints from the micropython input objects.
-    int a = mp_obj_get_int(a_obj);
-    int b = mp_obj_get_int(b_obj);
-
-    // Calculate the addition and convert to MicroPython object.
-    return mp_obj_new_int(a + b);
-}
-
-
-static mp_obj_t eigenmath_cmd() {
-	run_stdin();
-    // Calculate the addition and convert to MicroPython object.
-    return mp_obj_new_int(0);
-}
-
-static MP_DEFINE_CONST_FUN_OBJ_0(eigenmath_cmd_obj, eigenmath_cmd);
-
-// Define a Python reference to the function above.
-static MP_DEFINE_CONST_FUN_OBJ_2(example_add_ints_obj, example_add_ints);
-
-// This structure represents Timer instance objects.
-typedef struct _example_Timer_obj_t {
-    // All objects start with the base.
-    mp_obj_base_t base;
-    // Everything below can be thought of as instance attributes, but they
-    // cannot be accessed by MicroPython code directly. In this example we
-    // store the time at which the object was created.
-    mp_uint_t start_time;
-} example_Timer_obj_t;
-
-// This is the Timer.time() method. After creating a Timer object, this
-// can be called to get the time elapsed since creating the Timer.
-static mp_obj_t example_Timer_time(mp_obj_t self_in) {
-    // The first argument is self. It is cast to the *example_Timer_obj_t
-    // type so we can read its attributes.
-    example_Timer_obj_t *self = MP_OBJ_TO_PTR(self_in);
-
-    // Get the elapsed time and return it as a MicroPython integer.
-    mp_uint_t elapsed = mp_hal_ticks_ms() - self->start_time;
-    return mp_obj_new_int_from_uint(elapsed);
-}
-static MP_DEFINE_CONST_FUN_OBJ_1(example_Timer_time_obj, example_Timer_time);
-
-// This represents Timer.__new__ and Timer.__init__, which is called when
-// the user instantiates a Timer object.
-static mp_obj_t example_Timer_make_new(const mp_obj_type_t *type, size_t n_args, size_t n_kw, const mp_obj_t *args) {
-    // Allocates the new object and sets the type.
-    example_Timer_obj_t *self = mp_obj_malloc(example_Timer_obj_t, type);
-
-    // Initializes the time for this Timer instance.
-    self->start_time = mp_hal_ticks_ms();
-
-    // The make_new function always returns self.
-    return MP_OBJ_FROM_PTR(self);
-}
-
-// This collects all methods and other static class attributes of the Timer.
-// The table structure is similar to the module table, as detailed below.
-static const mp_rom_map_elem_t example_Timer_locals_dict_table[] = {
-    { MP_ROM_QSTR(MP_QSTR_time), MP_ROM_PTR(&example_Timer_time_obj) },
-};
-static MP_DEFINE_CONST_DICT(example_Timer_locals_dict, example_Timer_locals_dict_table);
-
-// This defines the type(Timer) object.
-MP_DEFINE_CONST_OBJ_TYPE(
-    example_type_Timer,
-    MP_QSTR_Timer,
-    MP_TYPE_FLAG_NONE,
-    make_new, example_Timer_make_new,
-    locals_dict, &example_Timer_locals_dict
-    );
-
-// What follows is a *separate* class definition that demonstrates more
-// advanced techniques to implement other Python-like features, such as:
-//
-// - A custom representation for __repr__ and __str__.
-// - Custom attribute handling to create a read/write "property".
-//
-// It reuses some of the elements of the basic Timer class. This is allowed
-// because they both use example_Timer_obj_t as the instance structure.
-
-// Handles AdvancedTimer.__repr__, AdvancedTimer.__str__.
-static void example_AdvancedTimer_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t kind) {
-
-    // Get the elapsed time. In this case, it's also a demonstration of calling
-    // the equivalent of self.time() in the C API. This is not usually very
-    // efficient, but it can sometimes be useful.
-    mp_uint_t elapsed = mp_obj_get_int(example_Timer_time(self_in));
-
-    // We'll make all representations print at least the class name.
-    mp_printf(print, "%q()", (qstr)MP_QSTR_AdvancedTimer);
-
-    // Decide what else to print based on print kind.
-    if (kind == PRINT_STR) {
-        // For __str__, let's attempt to make it more readable.
-        mp_printf(print, "  # created %d seconds ago", (int)(elapsed / 1000));
-    }
-}
-
-// Handles AdvancedTimer.seconds for reading and writing.
-static void example_AdvancedTimer_attribute_handler(mp_obj_t self_in, qstr attr, mp_obj_t *dest) {
-
-    // In this example, we only want to handle the .seconds attribute in a
-    // special way.
-    if (attr != MP_QSTR_seconds) {
-        // Attribute not found, continue lookup in locals dict. This way,
-        // methods like .time() will be handled normally.
-        dest[1] = MP_OBJ_SENTINEL;
-        return;
-    }
-
-    // Get reference to AdvancedTimer instance.
-    example_Timer_obj_t *self = MP_OBJ_TO_PTR(self_in);
-
-    // Check if this is a read operation.
-    if (dest[0] == MP_OBJ_NULL) {
-        // It's read, so "return" elapsed seconds by storing it in dest[0].
-        mp_uint_t elapsed = mp_hal_ticks_ms() - self->start_time;
-        dest[0] = mp_obj_new_int_from_uint(elapsed / 1000);
-        return;
-    }
-    // Check if this is a delete or store operation.
-    else if (dest[0] == MP_OBJ_SENTINEL) {
-        // It's delete or store. Now check which one.
-        if (dest[1] == MP_OBJ_NULL) {
-            // It's delete. But in this example we don't want to allow it
-            // so we just return.
-            return;
-        } else {
-            // It's write. First, get the value that the user is trying to set.
-            mp_uint_t desired_ms = mp_obj_get_int(dest[1]) * 1000;
-            // Use it to update the start time. This way, the next read will
-            // report the updated time.
-            self->start_time = mp_hal_ticks_ms() - desired_ms;
-
-            // Indicate successful store.
-            dest[0] = MP_OBJ_NULL;
-            return;
-        }
-    }
-}
-
-// This defines the type(AdvancedTimer) object.
-MP_DEFINE_CONST_OBJ_TYPE(
-    example_type_AdvancedTimer,
-    MP_QSTR_AdvancedTimer,
-    MP_TYPE_FLAG_NONE,
-    attr, example_AdvancedTimer_attribute_handler,
-    print, example_AdvancedTimer_print,
-    make_new, example_Timer_make_new,
-    locals_dict, &example_Timer_locals_dict
-    );
-
-// Define all attributes of the module.
-// Table entries are key/value pairs of the attribute name (a string)
-// and the MicroPython object reference.
-// All identifiers and strings are written as MP_QSTR_xxx and will be
-// optimized to word-sized integers by the build system (interned strings).
-static const mp_rom_map_elem_t eigenmath_module_globals_table[] = {
-    { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_eigenmath) },
-    { MP_ROM_QSTR(MP_QSTR_add_ints), MP_ROM_PTR(&example_add_ints_obj) },
-	{ MP_ROM_QSTR(MP_QSTR_cmd), MP_ROM_PTR(&eigenmath_cmd_obj) },
-    { MP_ROM_QSTR(MP_QSTR_Timer),    MP_ROM_PTR(&example_type_Timer) },
-    { MP_ROM_QSTR(MP_QSTR_AdvancedTimer),    MP_ROM_PTR(&example_type_AdvancedTimer) },
-};
-static MP_DEFINE_CONST_DICT(eigenmath_module_globals, eigenmath_module_globals_table);
-
-// Define module object.
-const mp_obj_module_t eigenmath_user_cmodule = {
-    .base = { &mp_type_module },
-    .globals = (mp_obj_dict_t *)&eigenmath_module_globals,
-};
-
-// Register the module to make it available in Python.
-MP_REGISTER_MODULE(MP_QSTR_eigenmath, eigenmath_user_cmodule);
\ No newline at end of file
diff --git a/user_modules/eigenmath/eigenmath.h b/user_modules/eigenmath/eigenmath.h
new file mode 100644
index 000000000..63e6b1698
--- /dev/null
+++ b/user_modules/eigenmath/eigenmath.h
@@ -0,0 +1,35 @@
+
+#ifndef EIGENMATH_H
+#define EIGENMATH_H
+
+#include <stddef.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+#define STRBUFLEN 1000
+#define BUCKETSIZE 100
+#define MAXDIM 24
+
+//extern struct atom *mem ;
+//extern struct atom **stack ; // 
+//extern struct atom **symtab ; // symbol table
+//extern struct atom **binding ;
+//extern struct atom **usrfunc ;
+//extern char *strbuf ;
+
+//extern uint32_t STACKSIZE ; // evaluation stack
+//extern uint32_t MAXATOMS ; // 10,240 atoms
+
+extern bool noprint;
+extern char *outbuf;
+extern int outbuf_index;
+extern void eigenmath_init(uint8_t *pHeap,size_t heapSize);
+extern void run(char *buf);
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* EHEAP_H */
\ No newline at end of file
diff --git a/user_modules/eigenmath/eigenmath_mpy.c b/user_modules/eigenmath/eigenmath_mpy.c
new file mode 100644
index 000000000..9d75426ec
--- /dev/null
+++ b/user_modules/eigenmath/eigenmath_mpy.c
@@ -0,0 +1,281 @@
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <ctype.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <string.h>
+#include <setjmp.h>
+#include <math.h>
+#include <errno.h>
+#include "py/obj.h"
+//#include "py/mpconfig.h"
+#include "py/misc.h"
+#include "py/runtime.h"
+#include "py/objstr.h"
+#include "shared/readline/readline.h"
+#include "py/binary.h"
+#include "py/gc.h"
+#include "py/stream.h"
+#include "eigenmath.h"
+#include "eheap.h"
+//-DPICO_STACK_SIZE=0x4000 ??
+
+typedef struct _mp_obj_eigenmath_t {
+    mp_obj_base_t base;
+    size_t heapSize;
+    uint8_t  *pHeap;
+} mp_obj_eigenmath_t;
+
+static void eigenmath_print(const mp_print_t *print, mp_obj_t self_in, mp_print_kind_t kind) {
+    mp_printf(print, "<EigenMath instance>");
+}
+
+static mp_obj_t eigenmath_make_new(const mp_obj_type_t *type,
+	size_t n_args, size_t n_kw,
+	const mp_obj_t *args) {
+
+	mp_arg_check_num(n_args, n_kw, 1, 1, false);
+	mp_obj_eigenmath_t *self = mp_obj_malloc(mp_obj_eigenmath_t, type);
+	self->base.type = type;
+    self->heapSize = mp_obj_get_int(args[0]); // 350 * 1024; // 350KB
+    //mp_printf(&mp_plat_print,"heapSize = %d\n", self->heapSize);
+
+    self->pHeap = (uint8_t *)m_malloc(self->heapSize);
+
+    //mp_printf(&mp_plat_print,"ptemp = %x\n", (uint32_t)(ptemp));
+    //mp_printf(&mp_plat_print,"self->pHeap = %x\n", (uint32_t)(self->pHeap));
+	if (self->pHeap == NULL){
+		mp_raise_msg(&mp_type_MemoryError, MP_ERROR_TEXT("Failed to initialize heap"));
+		return MP_OBJ_NULL;
+	} 
+    
+
+    eigenmath_init(self->pHeap,self->heapSize);
+
+	return MP_OBJ_FROM_PTR(self);
+}
+
+static mp_obj_t eigenmath_run(size_t n_args, const mp_obj_t *args) {
+	//mp_obj_eigenmath_t *self = MP_OBJ_TO_PTR(self_in);mp_obj_t input_str_obj
+	size_t len;
+    if (n_args >= 3){
+    mp_obj_t arg = args[2];
+    if (mp_obj_is_bool(arg) || mp_obj_is_int(arg)) {
+        noprint = mp_obj_is_true(arg);
+        } else {
+            mp_raise_TypeError(MP_ERROR_TEXT("expected a bool"));
+        }
+    } else {
+        noprint = false; 
+    }
+
+    if (!mp_obj_is_str(args[1])) {
+        mp_raise_TypeError(MP_ERROR_TEXT("expected a string as input"));
+    }
+    const char *buf = mp_obj_str_get_data(args[1], &len);
+
+	//GET_STR_DATA_LEN(input_str_obj, str, str_len);
+	run((char *)buf);
+
+    if (noprint == true){
+        return mp_obj_new_bytearray_by_ref(outbuf_index-1, outbuf);
+        // return memoryview
+        //return mp_obj_new_memoryview(BYTEARRAY_TYPECODE, outbuf);
+
+    }else{
+        return mp_const_none;
+    }
+	
+
+}
+static MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(eigenmath_run_obj,2,3, eigenmath_run);
+
+
+
+
+static mp_obj_t eigenmath_call(mp_obj_t self_in, size_t n_args, size_t n_kw, const mp_obj_t *args) {
+    //mp_obj_eigenmath_t *self = MP_OBJ_TO_PTR(self_in);
+    noprint = false;
+    if (n_args != 1 || n_kw != 0) {
+        mp_raise_TypeError(MP_ERROR_TEXT("Expected 1 positional argument"));
+    }
+
+    const char *cmd = mp_obj_str_get_str(args[0]);
+    run((char *)cmd);  // 
+    return mp_const_none;
+}
+
+static mp_obj_t eigenmath_cmd(mp_obj_t self_in) {
+	//mp_obj_eigenmath_t *self = MP_OBJ_TO_PTR(self_in);
+    vstr_t* vstr_inbuf = vstr_new(1);
+    for (;;) {
+        vstr_clear(vstr_inbuf);
+		int res = readline(vstr_inbuf,"eigenmath> ");
+		mp_printf(&mp_plat_print, "Eigenmath run:\n");
+		mp_printf(&mp_plat_print, "res=%d\n", res);
+		mp_printf(&mp_plat_print, "%s\n", vstr_inbuf->buf);
+		run(vstr_inbuf->buf);
+	}
+	return mp_const_none;
+
+}
+static MP_DEFINE_CONST_FUN_OBJ_1(eigenmath_cmd_obj, eigenmath_cmd);
+
+
+
+static mp_obj_t eigenmath_runfile(size_t n_args, const mp_obj_t *args ) {//mp_obj_t input_file_obj
+
+    if (n_args >= 3){
+        mp_obj_t arg = args[2];
+        if (mp_obj_is_bool(arg) || mp_obj_is_int(arg)) {
+            noprint = mp_obj_is_true(arg);
+        } else {
+            mp_raise_TypeError(MP_ERROR_TEXT("expected a bool"));
+        }
+    } else {
+         noprint = false; 
+    }
+    
+
+    const mp_stream_p_t *stream_p = mp_get_stream_raise(args[1], MP_STREAM_OP_READ | MP_STREAM_OP_IOCTL);
+    if (stream_p == NULL) {
+        mp_raise_TypeError(MP_ERROR_TEXT("expected a file-like object"));
+    }
+    int error = 0;
+
+    // get file size
+    struct mp_stream_seek_t seek = {
+        .offset = 0,
+        .whence = MP_SEEK_END,
+    };
+    mp_obj_t input_file_obj = args[1];
+    if (stream_p->ioctl(input_file_obj, MP_STREAM_SEEK, (uintptr_t)&seek, &error) == MP_STREAM_ERROR) {
+        mp_raise_OSError(error);
+    }
+    mp_off_t size = seek.offset;
+
+    // move to front
+    seek.offset = 0;
+    seek.whence = MP_SEEK_SET;
+    if (stream_p->ioctl(input_file_obj, MP_STREAM_SEEK, (uintptr_t)&seek, &error) == MP_STREAM_ERROR) {
+        mp_raise_OSError(error);
+    }
+
+    // get buffer
+    char *buf = m_new(char, size + 1);
+
+    // read file
+    mp_uint_t out_sz = stream_p->read(input_file_obj, buf, size, &error);
+    if (error != 0 || out_sz != size) {
+        m_del(char, buf, size + 1);
+        mp_raise_OSError(error);
+    }
+
+    // add end
+    buf[out_sz] = '\0';
+
+    // run
+    run(buf);
+
+    // release buffer
+    m_del(char, buf, size + 1);
+
+    if (noprint == true){
+        return mp_obj_new_bytearray_by_ref(outbuf_index-1, outbuf);
+        // return memoryview
+        //return mp_obj_new_memoryview(BYTEARRAY_TYPECODE, bytearray);
+
+    }else{
+        return mp_const_none;
+    }
+}
+static MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(eigenmath_runfile_obj, 2,3,eigenmath_runfile);
+
+
+
+
+
+extern int free_count;
+extern int MAXATOMS;
+static mp_obj_t eigenmath_status(mp_obj_t self_in) {
+	//mp_obj_eigenmath_t *self = MP_OBJ_TO_PTR(self_in);
+	int fragmentation = e_heap_fragmentation();
+    size_t free_bytes = e_heap_free();
+    size_t min_free = e_heap_min_free();
+    int num_atoms = free_count;
+    mp_printf(&mp_plat_print,"Heap fragmentation: %d%%\n", fragmentation);
+    mp_printf(&mp_plat_print,"Free bytes in Heap: %d\n", (int)free_bytes);
+    mp_printf(&mp_plat_print,"Minimum free bytes in Heap: %d\n", (int)min_free);
+    mp_printf(&mp_plat_print,"Number of free atoms: %d of %d\n", num_atoms,MAXATOMS);
+	return mp_const_none;
+
+}
+static MP_DEFINE_CONST_FUN_OBJ_1(eigenmath_status_obj, eigenmath_status);
+
+
+
+extern struct atom *zero;
+static mp_obj_t eigenmath_del(mp_obj_t self_in) {
+	mp_obj_eigenmath_t *self = MP_OBJ_TO_PTR(self_in);
+	m_free(&self->pHeap); // deinitialize the hea
+	zero = NULL;
+	return mp_const_none;
+}
+static MP_DEFINE_CONST_FUN_OBJ_1(eigenmath_del_obj, eigenmath_del);
+
+extern struct atom *zero;
+static mp_obj_t eigenmath_reset(mp_obj_t self_in) {
+    mp_obj_eigenmath_t *self = MP_OBJ_TO_PTR(self_in);
+	eigenmath_init(self->pHeap,self->heapSize);
+    zero = NULL;//triger the symbol table initialization
+    return mp_const_none;
+}
+static MP_DEFINE_CONST_FUN_OBJ_1(eigenmath_reset_obj, eigenmath_reset);
+
+mp_obj_t eigenmath_attr(mp_obj_t self_in, qstr attr, mp_obj_t *dest) {
+    if (dest[0] == MP_OBJ_NULL && attr == MP_QSTR___del__) {
+        dest[0] = MP_OBJ_FROM_PTR(&eigenmath_del_obj);
+        dest[1] = self_in;
+    }else{
+		// For any other attribute, indicate that lookup should continue in the locals dict
+		dest[1] = MP_OBJ_SENTINEL;
+		return MP_OBJ_NULL;
+		}
+
+    return MP_OBJ_NULL; 
+}
+static const mp_rom_map_elem_t eigenmath_locals_dict_table[] = {
+	{ MP_ROM_QSTR(MP_QSTR_run), MP_ROM_PTR(&eigenmath_run_obj) },
+    { MP_ROM_QSTR(MP_QSTR_runfile), MP_ROM_PTR(&eigenmath_runfile_obj) },
+    { MP_ROM_QSTR(MP_QSTR_cmd), MP_ROM_PTR(&eigenmath_cmd_obj) },
+	{ MP_ROM_QSTR(MP_QSTR_reset), MP_ROM_PTR(&eigenmath_reset_obj) },
+    { MP_ROM_QSTR(MP_QSTR_status), MP_ROM_PTR(&eigenmath_status_obj) },
+};
+static MP_DEFINE_CONST_DICT(eigenmath_locals_dict, eigenmath_locals_dict_table);
+
+
+MP_DEFINE_CONST_OBJ_TYPE(
+    eigenmath_type,
+    MP_QSTR_EigenMath,
+    MP_TYPE_FLAG_NONE,
+    make_new, eigenmath_make_new,
+    call,eigenmath_call, // call handler for the run method
+    attr, eigenmath_attr,           // attr handler before locals_dict
+    locals_dict, &eigenmath_locals_dict,
+    print, eigenmath_print
+);
+
+static const mp_rom_map_elem_t eigenmath_module_globals_table[] = {
+    { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_eigenmath) },
+    { MP_ROM_QSTR(MP_QSTR_EigenMath), MP_ROM_PTR(&eigenmath_type) },
+    };
+    static MP_DEFINE_CONST_DICT(mp_module_eigenmath_globals, eigenmath_module_globals_table);
+    
+    const mp_obj_module_t eigenmath_user_cmodule = {
+    .base = { &mp_type_module },
+    .globals = (mp_obj_dict_t*)&mp_module_eigenmath_globals,
+    };
+    
+    MP_REGISTER_MODULE(MP_QSTR_eigenmath, eigenmath_user_cmodule);
\ No newline at end of file
diff --git a/user_modules/eigenmath/micropython.cmake b/user_modules/eigenmath/micropython.cmake
index 0631eb107..b6819f75c 100644
--- a/user_modules/eigenmath/micropython.cmake
+++ b/user_modules/eigenmath/micropython.cmake
@@ -4,12 +4,16 @@ add_library(usermod_eigenmath INTERFACE)
 # Add our source files to the lib
 target_sources(usermod_eigenmath INTERFACE
     ${CMAKE_CURRENT_LIST_DIR}/eigenmath.c
+    ${CMAKE_CURRENT_LIST_DIR}/eheap.c
+    ${CMAKE_CURRENT_LIST_DIR}/eigenmath_mpy.c
 )
 
 # Add the current directory as an include directory.
 target_include_directories(usermod_eigenmath INTERFACE
     ${CMAKE_CURRENT_LIST_DIR}
 )
+set(PICO_STACK_SIZE 0x4000 CACHE STRING "App stack size" FORCE)
 
 # Link our INTERFACE library to the usermod target.
-target_link_libraries(usermod INTERFACE usermod_eigenmath)
\ No newline at end of file
+target_link_libraries(usermod INTERFACE usermod_eigenmath)
+
diff --git a/user_modules/eigenmath/micropython.mk b/user_modules/eigenmath/micropython.mk
index 577d7e1ba..dc5f9c5e9 100644
--- a/user_modules/eigenmath/micropython.mk
+++ b/user_modules/eigenmath/micropython.mk
@@ -1,8 +1,10 @@
-CEXAMPLE_MOD_DIR := $(USERMOD_DIR)
+PICOCALCDISPLAY_MOD_DIR := $(USERMOD_DIR)
 
 # Add all C files to SRC_USERMOD.
-SRC_USERMOD += $(CEXAMPLE_MOD_DIR)/eigenmath.c
+SRC_USERMOD += $(PICOCALCDISPLAY_MOD_DIR)/eigenmath.c\
+			  $(PICOCALCDISPLAY_MOD_DIR)/eigenmath_mpy.h\
+			  $(PICOCALCDISPLAY_MOD_DIR)/eheap.c\
 
 # We can add our module folder to include paths if needed
 # This is not actually needed in this example.
-CFLAGS_USERMOD += -I$(CEXAMPLE_MOD_DIR)
\ No newline at end of file
+CFLAGS_USERMOD += -I$(PICOCALCDISPLAY_MOD_DIR)