A ~10x speedup in the solver. Simplify equations before evaluating

or taking partials (constant folding). Also keep a little hash
table to mark with params are used in each equation, in order to
quickly discard trivial partial derivatives. This is solving a
64x64 system in <20 ms. I suspect this is now much faster than
Sketchflat.

Slightly fake situation, though, since substitution solver has not
yet been written, and no partitioning. I'll do those next.

[git-p4: depot-paths = "//depot/solvespace/": change = 1698]

expr.cpp

@@ -186,7 +186,7 @@ Expr *Expr::PartialWrt(hParam p) {
     Expr *da, *db;
 
     switch(op) {
-        case PARAM_PTR: oops();
+        case PARAM_PTR: return FromConstant(p.v == x.parp->h.v ? 1 : 0);
         case PARAM:     return FromConstant(p.v == x.parh.v ? 1 : 0);
 
         case CONSTANT:  return FromConstant(0);
@@ -218,6 +218,86 @@ Expr *Expr::PartialWrt(hParam p) {
     }
 }
 
+DWORD Expr::ParamsUsed(void) {
+    DWORD r = 0;
+    if(op == PARAM)     r |= (1 << (x.parh.v % 31));
+    if(op == PARAM_PTR) r |= (1 << (x.parp->h.v % 31));
+
+    int c = Children();
+    if(c >= 1)          r |= a->ParamsUsed();
+    if(c >= 2)          r |= b->ParamsUsed();
+    return r;
+}
+
+bool Expr::Tol(double a, double b) {
+    return fabs(a - b) < 0.001;
+}
+Expr *Expr::FoldConstants(void) {
+    Expr *n = AllocExpr();
+    *n = *this;
+
+    int c = Children();
+    if(c >= 1) n->a = a->FoldConstants();
+    if(c >= 2) n->b = b->FoldConstants();
+
+    switch(op) {
+        case PARAM_PTR:
+        case PARAM:
+        case CONSTANT:
+            break;
+
+        case MINUS:
+        case TIMES:
+        case DIV:
+        case PLUS:
+            // If both ops are known, then we can evaluate immediately
+            if(n->a->op == CONSTANT && n->b->op == CONSTANT) {
+                double nv = n->Eval();
+                n->op = CONSTANT;
+                n->x.v = nv;
+                break;
+            }
+            // x + 0 = 0 + x = x
+            if(op == PLUS && n->b->op == CONSTANT && Tol(n->b->x.v, 0)) {
+                *n = *(n->a); break;
+            }
+            if(op == PLUS && n->a->op == CONSTANT && Tol(n->a->x.v, 0)) {
+                *n = *(n->b); break;
+            }
+            // 1*x = x*1 = x
+            if(op == TIMES && n->b->op == CONSTANT && Tol(n->b->x.v, 1)) {
+                *n = *(n->a); break;
+            }
+            if(op == TIMES && n->a->op == CONSTANT && Tol(n->a->x.v, 1)) {
+                *n = *(n->b); break;
+            }
+            // 0*x = x*0 = 0
+            if(op == TIMES && n->b->op == CONSTANT && Tol(n->b->x.v, 0)) {
+                n->op = CONSTANT; n->x.v = 0; break;
+            }
+            if(op == TIMES && n->a->op == CONSTANT && Tol(n->a->x.v, 0)) {
+                n->op = CONSTANT; n->x.v = 0; break;
+            }
+
+            break;
+
+        case SQRT:
+        case SQUARE:
+        case NEGATE:
+        case SIN:
+        case COS:
+            if(n->a->op == CONSTANT) {
+                double nv = n->Eval();
+                n->op = CONSTANT;
+                n->x.v = nv;
+            }
+            break;
+
+        default: oops();
+    }
+    return n;
+}
+
 static char StringBuffer[4096];
 void Expr::App(char *s, ...) {
     va_list f;

expr.h

@@ -73,6 +73,9 @@ public:
 
     Expr *PartialWrt(hParam p);
     double Eval(void);
+    DWORD ParamsUsed(void);
+    static bool Tol(double a, double b);
+    Expr *FoldConstants(void);
 
     void ParamsToPointers(void);
 

file.cpp

@@ -11,6 +11,7 @@ void SolveSpace::NewFile(void) {
     // Our initial group, that contains the references.
     Group g;
     memset(&g, 0, sizeof(g));
+    g.visible = true;
     g.name.strcpy("#references");
     g.type = Group::DRAWING;
     g.h = Group::HGROUP_REFERENCES;

system.cpp

@@ -18,12 +18,22 @@ void System::WriteJacobian(int eqTag, int paramTag) {
         if(e->tag != eqTag) continue;
 
         mat.eq[i] = e->h;
-        mat.B.sym[i] = e->e->DeepCopyWithParamsAsPointers(&param, &(SS.param));
+        Expr *f = e->e->DeepCopyWithParamsAsPointers(&param, &(SS.param));
+        f = f->FoldConstants();
+
+        // Hash table (31 bits) to accelerate generation of zero partials.
+        DWORD scoreboard = f->ParamsUsed();
         for(j = 0; j < mat.n; j++) {
-            Expr *pd = e->e->PartialWrt(mat.param[j]);
-            mat.A.sym[i][j] = 
-                pd->DeepCopyWithParamsAsPointers(&param, &(SS.param));
+            Expr *pd;
+            if(scoreboard & (1 << (mat.param[j].v % 31))) { 
+                pd = f->PartialWrt(mat.param[j]);
+                pd = pd->FoldConstants();
+            } else {
+                pd = Expr::FromConstant(0);
+            }
+            mat.A.sym[i][j] = pd;
         }
+        mat.B.sym[i] = f;
         i++;
     }
     mat.m = i;
@@ -207,7 +217,7 @@ bool System::NewtonSolve(int tag) {
 
 bool System::Solve(void) {
     int i, j;
-
+    
 /*
     dbp("%d equations", eq.n);
     for(i = 0; i < eq.n; i++) {
@@ -220,7 +230,12 @@ bool System::Solve(void) {
     
     WriteJacobian(0, 0);
     EvalJacobian();
-/*
+
+/*    dbp("write/eval jacboian=%d", GetMilliseconds() - in);
+    for(i = 0; i < mat.m; i++) {
+        dbp("function %d: %s", i, mat.B.sym[i]->Print());
+    }
+    dbp("m=%d", mat.m);
     for(i = 0; i < mat.m; i++) {
         for(j = 0; j < mat.n; j++) {
             dbp("A[%d][%d] = %.3f", i, j, mat.A.num[i][j]);
@@ -245,7 +260,7 @@ bool System::Solve(void) {
         }
         p->tag = ASSUMED;
     }
-
+    
     bool ok = NewtonSolve(0);
 
     if(ok) {