diff --git a/expr.cpp b/expr.cpp
index cf75058..521f450 100644
--- a/expr.cpp
+++ b/expr.cpp
@@ -364,10 +364,10 @@ 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));
+QWORD Expr::ParamsUsed(void) {
+    QWORD r = 0;
+    if(op == PARAM)     r |= ((QWORD)1 << (x.parh.v % 61));
+    if(op == PARAM_PTR) r |= ((QWORD)1 << (x.parp->h.v % 61));
 
     int c = Children();
     if(c >= 1)          r |= a->ParamsUsed();
@@ -375,6 +375,16 @@ DWORD Expr::ParamsUsed(void) {
     return r;
 }
 
+bool Expr::DependsOn(hParam p) {
+    if(op == PARAM)     return (x.parh.v    == p.v);
+    if(op == PARAM_PTR) return (x.parp->h.v == p.v);
+
+    int c = Children();
+    if(c == 1)          return a->DependsOn(p);
+    if(c == 2)          return a->DependsOn(p) || b->DependsOn(p);
+    return false;
+}
+
 bool Expr::Tol(double a, double b) {
     return fabs(a - b) < 0.001;
 }
diff --git a/expr.h b/expr.h
index da0babe..06f47b2 100644
--- a/expr.h
+++ b/expr.h
@@ -77,7 +77,8 @@ public:
 
     Expr *PartialWrt(hParam p);
     double Eval(void);
-    DWORD ParamsUsed(void);
+    QWORD ParamsUsed(void);
+    bool DependsOn(hParam p);
     static bool Tol(double a, double b);
     Expr *FoldConstants(void);
     void Substitute(hParam oldh, hParam newh);
diff --git a/solvespace.h b/solvespace.h
index 060139d..af72bad 100644
--- a/solvespace.h
+++ b/solvespace.h
@@ -43,6 +43,7 @@ inline double WRAP_SYMMETRIC(double v, double n) {
 
 #define isforname(c) (isalnum(c) || (c) == '_' || (c) == '-' || (c) == '#')
 
+typedef unsigned __int64 QWORD;
 typedef signed long SDWORD;
 typedef signed short SWORD;
 
@@ -201,7 +202,7 @@ void MakePathAbsolute(char *base, char *path);
 
 class System {
 public:
-#define MAX_UNKNOWNS    200
+#define MAX_UNKNOWNS    1000
 
     EntityList                      entity;
     ParamList                       param;
diff --git a/system.cpp b/system.cpp
index 006ecb8..d026abf 100644
--- a/system.cpp
+++ b/system.cpp
@@ -24,11 +24,13 @@ void System::WriteJacobian(int tag) {
         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();
+        // Hash table (61 bits) to accelerate generation of zero partials.
+        QWORD scoreboard = f->ParamsUsed();
         for(j = 0; j < mat.n; j++) {
             Expr *pd;
-            if(scoreboard & (1 << (mat.param[j].v % 31))) { 
+            if(scoreboard & ((QWORD)1 << (mat.param[j].v % 61)) &&
+                f->DependsOn(mat.param[j]))
+            {
                 pd = f->PartialWrt(mat.param[j]);
                 pd = pd->FoldConstants();
                 pd = pd->DeepCopyWithParamsAsPointers(&param, &(SS.param));
@@ -225,7 +227,7 @@ bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS],
         for(ip = i+1; ip < n; ip++) {
             temp = A[ip][i]/A[i][i];
 
-            for(jp = 0; jp < n; jp++) {
+            for(jp = i; jp < n; jp++) {
                 A[ip][jp] -= temp*(A[i][jp]);
             }
             B[ip] -= temp*B[i];
@@ -291,7 +293,6 @@ bool System::SolveLeastSquares(void) {
 }
 
 bool System::NewtonSolve(int tag) {
-    WriteJacobian(tag);
     if(mat.m > mat.n) return false;
 
     int iter = 0;
@@ -438,6 +439,7 @@ void System::Solve(Group *g, bool andFindFree) {
 
         e->tag = alone;
         p->tag = alone;
+        WriteJacobian(alone);
         if(!NewtonSolve(alone)) {
             // Failed to converge, bail out early
             goto didnt_converge;