Speed up the inconsistent constraint detection; there's no reason

not to solve by substitution before rank testing. And report the
unsatisfied constraints when we don't converge.

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

sketch.h

@@ -554,6 +554,9 @@ public:
 class hEquation {
 public:
     DWORD v;
+
+    inline bool isFromConstraint(void);
+    inline hConstraint constraint(void);
 };
 
 class Equation {
@@ -599,5 +602,10 @@ inline hRequest hParam::request(void)
 inline hEquation hConstraint::equation(int i)
     { hEquation r; r.v = (v << 16) | i; return r; }
 
+inline bool hEquation::isFromConstraint(void)
+    { if(v & 0xc0000000) return false; else return true; }
+inline hConstraint hEquation::constraint(void)
+    { hConstraint r; r.v = (v >> 16); return r; }
+
 
 #endif

solvespace.h

@@ -212,7 +212,7 @@ public:
         }           B;
     } mat;
 
-    static const double RANK_MAG_TOLERANCE;
+    static const double RANK_MAG_TOLERANCE, CONVERGE_TOLERANCE;
     int CalculateRank(void);
     static bool SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS],
                                   double B[], int N);

system.cpp

@@ -1,6 +1,7 @@
 #include "solvespace.h"
 
 const double System::RANK_MAG_TOLERANCE = 1e-4;
+const double System::CONVERGE_TOLERANCE = 1e-10;
 
 void System::WriteJacobian(int tag) {
     int a, i, j;
@@ -328,7 +329,7 @@ bool System::NewtonSolve(int tag) {
             if(isnan(mat.B.num[i])) {
                 return false;
             }
-            if(fabs(mat.B.num[i]) > 1e-10) {
+            if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE) {
                 converged = false;
                 break;
             }
@@ -360,25 +361,39 @@ void System::WriteEquationsExceptFor(hConstraint hc, hGroup hg) {
 }
 
 void System::FindWhichToRemoveToFixJacobian(Group *g) {
-    int i;
+    int a, i;
     (g->solved.remove).Clear();
 
-    for(i = 0; i < SS.constraint.n; i++) {
+    for(a = 0; a < 2; a++) {
-        Constraint *c = &(SS.constraint.elem[i]);
+        for(i = 0; i < SS.constraint.n; i++) {
-        if(c->group.v != g->h.v) continue;
+            Constraint *c = &(SS.constraint.elem[i]);
+            if(c->group.v != g->h.v) continue;
+            if((c->type == Constraint::POINTS_COINCIDENT && a == 0) ||
+               (c->type != Constraint::POINTS_COINCIDENT && a == 1))
+            {
+                // Do the constraints in two passes: first everything but
+                // the point-coincident constraints, then only those
+                // constraints (so they appear last in the list).
+                continue;
+            }
 
-        param.ClearTags();
+            param.ClearTags();
-        eq.Clear();
+            eq.Clear();
-        WriteEquationsExceptFor(c->h, g->h);
+            WriteEquationsExceptFor(c->h, g->h);
-        eq.ClearTags();
+            eq.ClearTags();
 
-        WriteJacobian(0);
+            // It's a major speedup to solve the easy ones by substitution here,
-        EvalJacobian();
+            // and that doesn't break anything.
+            SolveBySubstitution();
 
-        int rank = CalculateRank();
+            WriteJacobian(0);
-        if(rank == mat.m) {
+            EvalJacobian();
-            // We fixed it by removing this constraint
+
-            (g->solved.remove).Add(&(c->h));
+            int rank = CalculateRank();
+            if(rank == mat.m) {
+                // We fixed it by removing this constraint
+                (g->solved.remove).Add(&(c->h));
+            }
         }
     }
 }
@@ -476,6 +491,25 @@ void System::Solve(Group *g) {
 
 didnt_converge:
     g->solved.how = Group::DIDNT_CONVERGE;
+    (g->solved.remove).Clear();
+    SS.constraint.ClearTags();
+    for(i = 0; i < eq.n; i++) {
+        if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE) {
+            // This constraint is unsatisfied.
+            if(!mat.eq[i].isFromConstraint()) continue;
+
+            hConstraint hc = mat.eq[i].constraint();
+            Constraint *c = SS.constraint.FindByIdNoOops(hc);
+            if(!c) continue;
+            // Don't double-show constraints that generated multiple
+            // unsatisfied equations
+            if(!c->tag) {
+                (g->solved.remove).Add(&(c->h));
+                c->tag = 1;
+            }
+        }
+    }
+
     TextWindow::ReportHowGroupSolved(g->h);
 }
 

textscreens.cpp

@@ -509,25 +509,26 @@ void TextWindow::ShowGroupSolveInfo(void) {
     Printf(true, "%FtGROUP   %E%s", g->DescriptionString());
     switch(g->solved.how) {
         case Group::DIDNT_CONVERGE:
-            Printf(true, "   %FxSOLVE FAILED!%Fd no convergence");
+            Printf(true, "%FxSOLVE FAILED!%Fd no convergence");
+            Printf(true, "the following constraints are unsatisfied");
             break;
 
-        case Group::SINGULAR_JACOBIAN: {
+        case Group::SINGULAR_JACOBIAN:
             Printf(true, "%FxSOLVE FAILED!%Fd inconsistent system");
             Printf(true, "remove any one of these to fix it");
-            for(int i = 0; i < g->solved.remove.n; i++) {
-                hConstraint hc = g->solved.remove.elem[i];
-                Constraint *c = SS.constraint.FindByIdNoOops(hc);
-                if(!c) continue;
-
-                Printf(false, "%Bp   %Fl%Ll%D%f%h%s%E",
-                    (i & 1) ? 'd' : 'a',
-                    c->h.v, (&TextWindow::ScreenSelectConstraint),
-                    (&TextWindow::ScreenHoverConstraint),
-                    c->DescriptionString());
-            }
             break;
-        }
+    }
+
+    for(int i = 0; i < g->solved.remove.n; i++) {
+        hConstraint hc = g->solved.remove.elem[i];
+        Constraint *c = SS.constraint.FindByIdNoOops(hc);
+        if(!c) continue;
+
+        Printf(false, "%Bp   %Fl%Ll%D%f%h%s%E",
+            (i & 1) ? 'd' : 'a',
+            c->h.v, (&TextWindow::ScreenSelectConstraint),
+            (&TextWindow::ScreenHoverConstraint),
+            c->DescriptionString());
     }
 
     Printf(true,  "It may be possible to fix the problem ");