Reimplement DivPivoting as DivProjected.

The old implementation was an approximation, whereas the new one is exact.

src/drawconstraint.cpp

@@ -148,7 +148,7 @@ int Constraint::DoLineTrimmedAgainstBox(Canvas *canvas, Canvas::hStroke hcs,
         }
         if(j < 4) continue;
 
-        double t = (p.Minus(a)).DivPivoting(dl);
+        double t = (p.Minus(a)).DivProjected(dl);
         tmin = min(t, tmin);
         tmax = max(t, tmax);
     }
@@ -642,7 +642,7 @@ void Constraint::DoLayout(DrawAs how, Canvas *canvas,
                 // Draw the projection marker from the closest point on the
                 // projected line to the projected point on the real line.
                 Vector lAB = (lA.Minus(lB));
-                double t = (lA.Minus(closest)).DivPivoting(lAB);
+                double t = (lA.Minus(closest)).DivProjected(lAB);
 
                 Vector lA = SK.GetEntity(line->point[0])->PointGetNum();
                 Vector lB = SK.GetEntity(line->point[1])->PointGetNum();

src/dsc.h

@@ -120,7 +120,7 @@ public:
     Vector ScaledBy(double s) const;
     Vector ProjectInto(hEntity wrkpl) const;
     Vector ProjectVectorInto(hEntity wrkpl) const;
-    double DivPivoting(Vector delta) const;
+    double DivProjected(Vector delta) const;
     Vector ClosestOrtho() const;
     void MakeMaxMin(Vector *maxv, Vector *minv) const;
     Vector ClampWithin(double minv, double maxv) const;
@@ -187,7 +187,7 @@ public:
     Point2d Plus(const Point2d &b) const;
     Point2d Minus(const Point2d &b) const;
     Point2d ScaledBy(double s) const;
-    double DivPivoting(Point2d delta) const;
+    double DivProjected(Point2d delta) const;
     double Dot(Point2d p) const;
     double DistanceTo(const Point2d &p) const;
     double DistanceToLine(const Point2d &p0, const Point2d &dp, bool asSegment) const;

src/modify.cpp

@@ -373,8 +373,8 @@ void GraphicsWindow::MakeTangentArc() {
         tp[1] = t[1];
 
         // And convert those points to parameter values along the curve.
-        t[0] += (pa0.Minus(p0)).DivPivoting(t0);
+        t[0] += (pa0.Minus(p0)).DivProjected(t0);
-        t[1] += (pa1.Minus(p1)).DivPivoting(t1);
+        t[1] += (pa1.Minus(p1)).DivProjected(t1);
     }
 
     // Stupid check for convergence, and for an out of range result (as

src/polygon.cpp

@@ -164,13 +164,13 @@ bool SEdge::EdgeCrosses(Vector ea, Vector eb, Vector *ppi, SPointList *spl) cons
         // on the other
         bool inters = false;
         double t;
-        t = a.Minus(ea).DivPivoting(d);
+        t = a.Minus(ea).DivProjected(d);
         if(t > t_eps && t < (1 - t_eps)) inters = true;
-        t = b.Minus(ea).DivPivoting(d);
+        t = b.Minus(ea).DivProjected(d);
         if(t > t_eps && t < (1 - t_eps)) inters = true;
-        t = ea.Minus(a).DivPivoting(dthis);
+        t = ea.Minus(a).DivProjected(dthis);
         if(t > tthis_eps && t < (1 - tthis_eps)) inters = true;
-        t = eb.Minus(a).DivPivoting(dthis);
+        t = eb.Minus(a).DivProjected(dthis);
         if(t > tthis_eps && t < (1 - tthis_eps)) inters = true;
 
         if(inters) {
@@ -500,8 +500,8 @@ void SEdgeList::MergeCollinearSegments(Vector a, Vector b) {
     const Vector lineStart = a;
     const Vector lineDirection = b.Minus(a);
     std::sort(l.begin(), l.end(), [&](const SEdge &a, const SEdge &b) {
-        double ta = (a.a.Minus(lineStart)).DivPivoting(lineDirection);
+        double ta = (a.a.Minus(lineStart)).DivProjected(lineDirection);
-        double tb = (b.a.Minus(lineStart)).DivPivoting(lineDirection);
+        double tb = (b.a.Minus(lineStart)).DivProjected(lineDirection);
 
         return (ta < tb);
     });

src/srf/boolean.cpp

@@ -96,8 +96,8 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
             const Vector lineStart     = prev.p;
             const Vector lineDirection = (p->p).Minus(prev.p);
             std::sort(il.begin(), il.end(), [&](const SInter &a, const SInter &b) {
-                double ta = (a.p.Minus(lineStart)).DivPivoting(lineDirection);
+                double ta = (a.p.Minus(lineStart)).DivProjected(lineDirection);
-                double tb = (b.p.Minus(lineStart)).DivPivoting(lineDirection);
+                double tb = (b.p.Minus(lineStart)).DivProjected(lineDirection);
 
                 return (ta < tb);
             });

src/srf/ratpoly.cpp

@@ -155,7 +155,7 @@ void SBezier::ClosestPointTo(Vector p, double *t, bool mustConverge) const {
 
         Vector dp = TangentAt(*t);
         Vector pc = p.ClosestPointOnLine(p0, dp);
-        *t += (pc.Minus(p0)).DivPivoting(dp);
+        *t += (pc.Minus(p0)).DivProjected(dp);
     }
     if(mustConverge) {
         dbp("didn't converge (closest point on bezier curve)");

src/srf/raycast.cpp

@@ -311,7 +311,7 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
         }
         int i;
         for(i = 0; i < ip_n; i++) {
-            double t = (ip[i].Minus(ap)).DivPivoting(bp.Minus(ap));
+            double t = (ip[i].Minus(ap)).DivProjected(bp.Minus(ap));
             // This is a point on the circle; but is it on the arc?
             Point2d pp = ap.Plus((bp.Minus(ap)).ScaledBy(t));
             double theta = atan2(pp.y, pp.x);
@@ -354,7 +354,7 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
 
         // Make sure the point lies within the finite line segment
         Vector pxyz = PointAt(puv.x, puv.y);
-        double t = (pxyz.Minus(a)).DivPivoting(ba);
+        double t = (pxyz.Minus(a)).DivProjected(ba);
         if(asSegment && (t > 1 - LENGTH_EPS/bam || t < LENGTH_EPS/bam)) {
             continue;
         }
@@ -566,7 +566,7 @@ bool SShell::ClassifyEdge(Class *indir, Class *outdir,
 
         SInter *si;
         for(si = l.First(); si; si = l.NextAfter(si)) {
-            double t = ((si->p).Minus(p)).DivPivoting(ray);
+            double t = ((si->p).Minus(p)).DivProjected(ray);
             if(t*ray.Magnitude() < -LENGTH_EPS) {
                 // wrong side, doesn't count
                 continue;

src/util.cpp

@@ -607,7 +607,7 @@ bool Vector::OnLineSegment(Vector a, Vector b, double tol) const {
 
     if(distsq >= tol*tol) return false;
 
-    double t = (this->Minus(a)).DivPivoting(d);
+    double t = (this->Minus(a)).DivProjected(d);
     // On-endpoint already tested
     if(t < 0 || t > 1) return false;
     return true;
@@ -696,16 +696,9 @@ Vector4 Vector::Project4d() const {
     return Vector4::From(1, x, y, z);
 }
 
-double Vector::DivPivoting(Vector delta) const {
+double Vector::DivProjected(Vector delta) const {
-    double mx = fabs(delta.x), my = fabs(delta.y), mz = fabs(delta.z);
+    return (x*delta.x + y*delta.y + z*delta.z)
-
+         / (delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
-    if(mx > my && mx > mz) {
-        return x/delta.x;
-    } else if(my > mz) {
-        return y/delta.y;
-    } else {
-        return z/delta.z;
-    }
 }
 
 Vector Vector::ClosestOrtho() const {
@@ -995,12 +988,8 @@ Point2d Point2d::ScaledBy(double s) const {
     return { x * s, y * s };
 }
 
-double Point2d::DivPivoting(Point2d delta) const {
+double Point2d::DivProjected(Point2d delta) const {
-    if(fabs(delta.x) > fabs(delta.y)) {
+    return (x*delta.x + y*delta.y) / (delta.x*delta.x + delta.y*delta.y);
-        return x/delta.x;
-    } else {
-        return y/delta.y;
-    }
 }
 
 double Point2d::MagSquared() const {