diff --git a/constraint.cpp b/constraint.cpp
index ee70440..181cdda 100644
--- a/constraint.cpp
+++ b/constraint.cpp
@@ -591,12 +591,12 @@ void Constraint::MenuConstrain(int id) {
                 }
                 Vector l0 = SK.GetEntity(line->point[0])->PointGetNum(),
                        l1 = SK.GetEntity(line->point[1])->PointGetNum();
-                Vector a0 = SK.GetEntity(cubic->point[0])->PointGetNum(),
-                       a3 = SK.GetEntity(cubic->point[3])->PointGetNum();
+                Vector as = cubic->CubicGetStartNum(),
+                       af = cubic->CubicGetFinishNum();
 
-                if(l0.Equals(a0) || l1.Equals(a0)) {
+                if(l0.Equals(as) || l1.Equals(as)) {
                     c.other = false;
-                } else if(l0.Equals(a3) || l1.Equals(a3)) {
+                } else if(l0.Equals(af) || l1.Equals(af)) {
                     c.other = true;
                 } else {
                     Error("The tangent cubic and line segment must share an "
diff --git a/constrainteq.cpp b/constrainteq.cpp
index e9d8c4a..fa6c3cc 100644
--- a/constrainteq.cpp
+++ b/constrainteq.cpp
@@ -639,12 +639,12 @@ void ConstraintBase::GenerateReal(IdList<Equation,hEquation> *l) {
             EntityBase *cubic = SK.GetEntity(entityA);
             EntityBase *line  = SK.GetEntity(entityB);
             
-            ExprVector endpoint =
-                SK.GetEntity(cubic->point[other ? 3 : 0])->PointGetExprs();
-            ExprVector ctrlpoint = 
-                SK.GetEntity(cubic->point[other ? 2 : 1])->PointGetExprs();
-            
-            ExprVector a = endpoint.Minus(ctrlpoint);
+            ExprVector a;
+            if(other) {
+                a = cubic->CubicGetFinishTangentExprs();
+            } else {
+                a = cubic->CubicGetStartTangentExprs();
+            }
 
             ExprVector b = line->VectorGetExprs();
 
diff --git a/draw.cpp b/draw.cpp
index 6992c9b..db18be5 100644
--- a/draw.cpp
+++ b/draw.cpp
@@ -172,6 +172,7 @@ void GraphicsWindow::GroupSelection(void) {
                 case Entity::WORKPLANE:     (gs.workplanes)++; break;
                 case Entity::LINE_SEGMENT:  (gs.lineSegments)++; break;
                 case Entity::CUBIC:         (gs.cubics)++; break;
+                // but don't count periodic cubics
 
                 case Entity::ARC_OF_CIRCLE:
                     (gs.circlesOrArcs)++;
@@ -203,7 +204,17 @@ void GraphicsWindow::HitTestMakeSelection(Point2d mp) {
     for(i = 0; i < SK.entity.n; i++) {
         Entity *e = &(SK.entity.elem[i]);
         // Don't hover whatever's being dragged.
-        if(e->h.request().v == pending.point.request().v) continue;
+        if(e->h.request().v == pending.point.request().v) {
+            // The one exception is when we're creating a new cubic; we
+            // want to be able to hover the first point, because that's
+            // how we turn it into a periodic spline.
+            if(!e->IsPoint()) continue;
+            if(!e->h.isFromRequest()) continue;
+            Request *r = SK.GetRequest(e->h.request());
+            if(r->type != Request::CUBIC) continue;
+            if(r->extraPoints < 2) continue;
+            if(e->h.v != r->h.entity(1).v) continue;
+        }
 
         d = e->GetDistance(mp);
         if(d < 10 && d < dmin) {
diff --git a/drawconstraint.cpp b/drawconstraint.cpp
index 9381c7f..b648288 100644
--- a/drawconstraint.cpp
+++ b/drawconstraint.cpp
@@ -699,8 +699,8 @@ void Constraint::DrawOrGetDistance(Vector *labelPos) {
                 }
 
                 Entity *cubic = SK.GetEntity(entityA);
-                Vector p =
-                    SK.GetEntity(cubic->point[other ? 3 : 0])->PointGetNum();
+                Vector p = other ? cubic->CubicGetFinishNum() :
+                                   cubic->CubicGetStartNum();  
                 Vector dir = SK.GetEntity(entityB)->VectorGetNum();
                 Vector out = n.Cross(dir);
                 textAt = p.Plus(out.WithMagnitude(14/SS.GW.scale));
diff --git a/drawentity.cpp b/drawentity.cpp
index 43a3766..c704358 100644
--- a/drawentity.cpp
+++ b/drawentity.cpp
@@ -215,6 +215,143 @@ bool Entity::PointIsFromReferences(void) {
     return h.request().IsFromReferences();
 }
 
+//-----------------------------------------------------------------------------
+// Compute a cubic, second derivative continuous, interpolating spline. Same
+// routine for periodic splines (in a loop) or open splines (with specified
+// end tangents).
+//-----------------------------------------------------------------------------
+void Entity::ComputeInterpolatingSpline(SBezierList *sbl, bool periodic) {
+    static const int MAX_N = BandedMatrix::MAX_UNKNOWNS;
+    int ep = extraPoints;
+
+    // The number of unknowns to solve for.
+    int n   = periodic ? 3 + ep : ep;
+    if(n >= MAX_N) oops();
+    // The number of on-curve points, one more than the number of segments.
+    int pts = periodic ? 4 + ep : 2 + ep;
+
+    int i, j, a;
+    
+    // The starting and finishing control points that define our end tangents
+    // (if the spline isn't periodic), and the on-curve points.
+    Vector ctrl_s, ctrl_f, pt[MAX_N+4];
+    if(periodic) {
+        for(i = 0; i < ep + 3; i++) {
+            pt[i] = SK.GetEntity(point[i])->PointGetNum();
+        }
+        pt[i++] = SK.GetEntity(point[0])->PointGetNum();
+    } else {
+        ctrl_s = SK.GetEntity(point[1])->PointGetNum();
+        ctrl_f = SK.GetEntity(point[ep+2])->PointGetNum();
+        j = 0;
+        pt[j++] = SK.GetEntity(point[0])->PointGetNum();
+        for(i = 2; i <= ep + 1; i++) {
+            pt[j++] = SK.GetEntity(point[i])->PointGetNum();
+        }
+        pt[j++] = SK.GetEntity(point[ep+3])->PointGetNum();
+    }
+
+    // The unknowns that we will be solving for, a set for each coordinate.
+    double Xx[MAX_N], Xy[MAX_N], Xz[MAX_N];
+    // For a cubic Bezier section f(t) as t goes from 0 to 1,
+    //    f' (0) = 3*(P1 - P0)
+    //    f' (1) = 3*(P3 - P2)
+    //    f''(0) = 6*(P0 - 2*P1 + P2)
+    //    f''(1) = 6*(P3 - 2*P2 + P1)
+    for(a = 0; a < 3; a++) {
+        BandedMatrix bm;
+        ZERO(&bm);
+        bm.n = n;
+   
+        for(i = 0; i < n; i++) {
+            int im, it, ip;
+            if(periodic) {
+                im = WRAP(i - 1, n);
+                it = i;
+                ip = WRAP(i + 1, n);
+            } else {
+                im = i;
+                it = i + 1;
+                ip = i + 2;
+            }
+            // All of these are expressed in terms of a constant part, and
+            // of X[i-1], X[i], and X[i+1]; so let these be the four
+            // components of that vector;
+            Vector4 A, B, C, D, E;
+            // The on-curve interpolated point
+            C = Vector4::From((pt[it]).Element(a), 0, 0, 0);
+            // control point one back, C - X[i]
+            B = C.Plus(Vector4::From(0, 0, -1, 0));
+            // control point one forward, C + X[i]
+            D = C.Plus(Vector4::From(0, 0, 1, 0));
+            // control point two back
+            if(i == 0 && !periodic) {
+                A = Vector4::From(ctrl_s.Element(a), 0, 0, 0);
+            } else {
+                // pt[im] + X[i-1]
+                A = Vector4::From(pt[im].Element(a), 1, 0, 0);
+            }
+            // control point two forward
+            if(i == (n - 1) && !periodic) {
+                E = Vector4::From(ctrl_f.Element(a), 0, 0, 0);
+            } else {
+                // pt[ip] - X[i+1]
+                E = Vector4::From((pt[ip]).Element(a), 0, 0, -1);
+            }
+            // Write the second derivatives of each segment, dropping constant
+            Vector4 fprev_pp = (C.Minus(B.ScaledBy(2))).Plus(A),
+                    fnext_pp = (C.Minus(D.ScaledBy(2))).Plus(E),
+                    eq       = fprev_pp.Minus(fnext_pp);
+
+            bm.B[i] = -eq.w;
+            if(periodic) {
+                bm.A[i][WRAP(i-2, n)] = eq.x;
+                bm.A[i][WRAP(i-1, n)] = eq.y;
+                bm.A[i][i]            = eq.z;
+            } else {
+                // The wrapping would work, except when n = 1 and everything
+                // wraps to zero...
+                if(i > 0)     bm.A[i][i - 1] = eq.x; 
+                              bm.A[i][i]     = eq.y;
+                if(i < (n-1)) bm.A[i][i + 1] = eq.z;
+            }
+        }
+        bm.Solve();
+        double *X = (a == 0) ? Xx :
+                    (a == 1) ? Xy :
+                               Xz;
+        memcpy(X, bm.X, n*sizeof(double));
+    }
+
+    for(i = 0; i < pts - 1; i++) {
+        Vector p0, p1, p2, p3;
+        if(periodic) {
+            p0 = pt[i];
+            int iw = WRAP(i - 1, n);
+            p1 = p0.Plus(Vector::From(Xx[iw], Xy[iw], Xz[iw]));
+        } else if(i == 0) {
+            p0 = pt[0];
+            p1 = ctrl_s;
+        } else {
+            p0 = pt[i];
+            p1 = p0.Plus(Vector::From(Xx[i-1], Xy[i-1], Xz[i-1]));
+        }
+        if(periodic) {
+            p3 = pt[i+1];
+            int iw = WRAP(i, n);
+            p2 = p3.Minus(Vector::From(Xx[iw], Xy[iw], Xz[iw]));
+        } else if(i == (pts - 2)) {
+            p3 = pt[pts-1];
+            p2 = ctrl_f;
+        } else {
+            p3 = pt[i+1];
+            p2 = p3.Minus(Vector::From(Xx[i], Xy[i], Xz[i]));
+        }
+        SBezier sb = SBezier::From(p0, p1, p2, p3);
+        sbl->l.Add(&sb);
+    }
+}
+
 void Entity::GenerateBezierCurves(SBezierList *sbl) {
     SBezier sb;
 
@@ -228,15 +365,13 @@ void Entity::GenerateBezierCurves(SBezierList *sbl) {
             sbl->l.Add(&sb);
             break;
         }
-        case CUBIC: {
-            Vector p0 = SK.GetEntity(point[0])->PointGetNum();
-            Vector p1 = SK.GetEntity(point[1])->PointGetNum();
-            Vector p2 = SK.GetEntity(point[2])->PointGetNum();
-            Vector p3 = SK.GetEntity(point[3])->PointGetNum();
-            sb = SBezier::From(p0, p1, p2, p3);
-            sbl->l.Add(&sb);
+        case CUBIC:
+            ComputeInterpolatingSpline(sbl, false);
+            break;
+
+        case CUBIC_PERIODIC:
+            ComputeInterpolatingSpline(sbl, true);
             break;
-        }
 
         case CIRCLE:
         case ARC_OF_CIRCLE: {
@@ -476,6 +611,7 @@ void Entity::DrawOrGetDistance(void) {
         case CIRCLE:
         case ARC_OF_CIRCLE:
         case CUBIC:
+        case CUBIC_PERIODIC:
         case TTF_TEXT:
             // Nothing but the curve(s).
             break;
diff --git a/dsc.h b/dsc.h
index 8e951d6..9c2727e 100644
--- a/dsc.h
+++ b/dsc.h
@@ -367,4 +367,18 @@ public:
     }
 };
 
+class BandedMatrix {
+public:
+    static const int MAX_UNKNOWNS   = 16;
+    static const int RIGHT_OF_DIAG  = 1;
+    static const int LEFT_OF_DIAG   = 2;
+
+    double A[MAX_UNKNOWNS][MAX_UNKNOWNS];
+    double B[MAX_UNKNOWNS];
+    double X[MAX_UNKNOWNS];
+    int n;
+
+    void Solve(void);
+};
+
 #endif
diff --git a/entity.cpp b/entity.cpp
index b2aeae2..4624e75 100644
--- a/entity.cpp
+++ b/entity.cpp
@@ -114,6 +114,23 @@ void EntityBase::ArcGetAngles(double *thetaa, double *thetab, double *dtheta) {
     while(*dtheta > (2*PI)) *dtheta -= 2*PI;
 }
 
+Vector EntityBase::CubicGetStartNum(void) {
+    return SK.GetEntity(point[0])->PointGetNum();
+}
+Vector EntityBase::CubicGetFinishNum(void) {
+    return SK.GetEntity(point[3+extraPoints])->PointGetNum();
+}
+ExprVector EntityBase::CubicGetStartTangentExprs(void) {
+    ExprVector pon  = SK.GetEntity(point[0])->PointGetExprs(),
+               poff = SK.GetEntity(point[1])->PointGetExprs();
+    return (pon.Minus(poff));
+}
+ExprVector EntityBase::CubicGetFinishTangentExprs(void) {
+    ExprVector pon  = SK.GetEntity(point[3+extraPoints])->PointGetExprs(),
+               poff = SK.GetEntity(point[2+extraPoints])->PointGetExprs();
+    return (pon.Minus(poff));
+}
+
 bool EntityBase::IsWorkplane(void) {
     return (type == WORKPLANE);
 }
diff --git a/file.cpp b/file.cpp
index f7e752f..12d774b 100644
--- a/file.cpp
+++ b/file.cpp
@@ -103,6 +103,7 @@ const SolveSpace::SaveTable SolveSpace::SAVED[] = {
 
     { 'r',  "Request.h.v",              'x',    &(SS.sv.r.h.v)                },
     { 'r',  "Request.type",             'd',    &(SS.sv.r.type)               },
+    { 'r',  "Request.extraPoints",      'd',    &(SS.sv.r.extraPoints)        },
     { 'r',  "Request.workplane.v",      'x',    &(SS.sv.r.workplane.v)        },
     { 'r',  "Request.group.v",          'x',    &(SS.sv.r.group.v)            },
     { 'r',  "Request.construction",     'b',    &(SS.sv.r.construction)       },
@@ -120,6 +121,15 @@ const SolveSpace::SaveTable SolveSpace::SAVED[] = {
     { 'e',  "Entity.point[1].v",        'x',    &(SS.sv.e.point[1].v)         },
     { 'e',  "Entity.point[2].v",        'x',    &(SS.sv.e.point[2].v)         },
     { 'e',  "Entity.point[3].v",        'x',    &(SS.sv.e.point[3].v)         },
+    { 'e',  "Entity.point[4].v",        'x',    &(SS.sv.e.point[4].v)         },
+    { 'e',  "Entity.point[5].v",        'x',    &(SS.sv.e.point[5].v)         },
+    { 'e',  "Entity.point[6].v",        'x',    &(SS.sv.e.point[6].v)         },
+    { 'e',  "Entity.point[7].v",        'x',    &(SS.sv.e.point[7].v)         },
+    { 'e',  "Entity.point[8].v",        'x',    &(SS.sv.e.point[8].v)         },
+    { 'e',  "Entity.point[9].v",        'x',    &(SS.sv.e.point[9].v)         },
+    { 'e',  "Entity.point[10].v",       'x',    &(SS.sv.e.point[10].v)        },
+    { 'e',  "Entity.point[11].v",       'x',    &(SS.sv.e.point[11].v)        },
+    { 'e',  "Entity.extraPoints",       'd',    &(SS.sv.e.extraPoints)        },
     { 'e',  "Entity.normal.v",          'x',    &(SS.sv.e.normal.v)           },
     { 'e',  "Entity.distance.v",        'x',    &(SS.sv.e.distance.v)         },
     { 'e',  "Entity.workplane.v",       'x',    &(SS.sv.e.workplane.v)        },
diff --git a/group.cpp b/group.cpp
index 972a5dd..62d200e 100644
--- a/group.cpp
+++ b/group.cpp
@@ -636,6 +636,7 @@ void Group::CopyEntity(IdList<Entity,hEntity> *el,
     Entity en;
     memset(&en, 0, sizeof(en));
     en.type = ep->type;
+    en.extraPoints = ep->extraPoints;
     en.h = Remap(ep->h, remap);
     en.timesApplied = timesApplied;
     en.group = h;
@@ -652,12 +653,20 @@ void Group::CopyEntity(IdList<Entity,hEntity> *el,
             en.point[1] = Remap(ep->point[1], remap);
             break;
 
-        case Entity::CUBIC:
-            en.point[0] = Remap(ep->point[0], remap);
-            en.point[1] = Remap(ep->point[1], remap);
-            en.point[2] = Remap(ep->point[2], remap);
-            en.point[3] = Remap(ep->point[3], remap);
+        case Entity::CUBIC: {
+            int i;
+            for(i = 0; i < 4 + ep->extraPoints; i++) {
+                en.point[i] = Remap(ep->point[i], remap);
+            }
             break;
+        }
+        case Entity::CUBIC_PERIODIC: {
+            int i;
+            for(i = 0; i < 3 + ep->extraPoints; i++) {
+                en.point[i] = Remap(ep->point[i], remap);
+            }
+            break;
+        }
 
         case Entity::CIRCLE:
             en.point[0] = Remap(ep->point[0], remap);
diff --git a/modify.cpp b/modify.cpp
index e5b0746..03963dc 100644
--- a/modify.cpp
+++ b/modify.cpp
@@ -310,7 +310,7 @@ hEntity GraphicsWindow::SplitCubic(hEntity he, Vector pinter) {
     b01.ClosestPointTo(pinter, &t, true);
     b01.SplitAt(t, &b0i, &bi1);
 
-    // Add the two line segments this one gets split into.
+    // Add the two cubic segments this one gets split into.
     hRequest r0i = AddRequest(Request::CUBIC, false),
              ri1 = AddRequest(Request::CUBIC, false);
     // Don't get entities till after adding, realloc issues
@@ -343,7 +343,7 @@ hEntity GraphicsWindow::SplitEntity(hEntity he, Vector pinter) {
         ret = SplitCircle(he, pinter);
     } else if(e->type == Entity::LINE_SEGMENT) {
         ret = SplitLine(he, pinter);
-    } else if(e->type == Entity::CUBIC) {
+    } else if(e->type == Entity::CUBIC && e->extraPoints == 0) {
         ret = SplitCubic(he, pinter);
     } else {
         Error("Couldn't split this entity; lines, circles, or cubics only.");
diff --git a/mouse.cpp b/mouse.cpp
index ace4162..cb768e9 100644
--- a/mouse.cpp
+++ b/mouse.cpp
@@ -205,12 +205,24 @@ void GraphicsWindow::MouseMoved(double x, double y, bool leftDown,
             HitTestMakeSelection(mp);
 
             hRequest hr = pending.point.request();
-            Vector p0 = SK.GetEntity(hr.entity(1))->PointGetNum();
-            Vector p3 = SK.GetEntity(hr.entity(4))->PointGetNum();
-            Vector p1 = p0.ScaledBy(2.0/3).Plus(p3.ScaledBy(1.0/3));
-            SK.GetEntity(hr.entity(2))->PointForceTo(p1);
-            Vector p2 = p0.ScaledBy(1.0/3).Plus(p3.ScaledBy(2.0/3));
-            SK.GetEntity(hr.entity(3))->PointForceTo(p2);
+            if(pending.point.v == hr.entity(4).v) {
+                // The very first segment; dragging final point drags both
+                // tangent points.
+                Vector p0 = SK.GetEntity(hr.entity(1))->PointGetNum(),
+                       p3 = SK.GetEntity(hr.entity(4))->PointGetNum(),
+                       p1 = p0.ScaledBy(2.0/3).Plus(p3.ScaledBy(1.0/3)),
+                       p2 = p0.ScaledBy(1.0/3).Plus(p3.ScaledBy(2.0/3));
+                SK.GetEntity(hr.entity(1+1))->PointForceTo(p1);
+                SK.GetEntity(hr.entity(1+2))->PointForceTo(p2);
+            } else {
+                // A subsequent segment; dragging point drags only final
+                // tangent point.
+                int i = SK.GetEntity(hr.entity(0))->extraPoints;
+                Vector pn   = SK.GetEntity(hr.entity(4+i))->PointGetNum(),
+                       pnm2 = SK.GetEntity(hr.entity(2+i))->PointGetNum(),
+                       pnm1 = (pn.Plus(pnm2)).ScaledBy(0.5);
+                SK.GetEntity(hr.entity(3+i))->PointForceTo(pnm1);
+            }
 
             SS.MarkGroupDirtyByEntity(pending.point);
             break;
@@ -282,6 +294,7 @@ void GraphicsWindow::MouseMoved(double x, double y, bool leftDown,
 
 void GraphicsWindow::ClearPending(void) {
     memset(&pending, 0, sizeof(pending));
+    SS.later.showTW = true;
 }
 
 void GraphicsWindow::MouseMiddleOrRightDown(double x, double y) {
@@ -319,15 +332,19 @@ void GraphicsWindow::MouseRightUp(double x, double y) {
     if(orig.startedMoving) return;
 
     if(context.active) return;
-    context.active = true;
 
-    if(pending.operation == DRAGGING_NEW_LINE_POINT) {
+    if(pending.operation == DRAGGING_NEW_LINE_POINT ||
+       pending.operation == DRAGGING_NEW_CUBIC_POINT)
+    {
         // Special case; use a right click to stop drawing lines, since
         // a left click would draw another one. This is quicker and more
-        // intuitive than hitting escape.
+        // intuitive than hitting escape. Likewise for new cubic segments.
         ClearPending();
+        return;
     }
 
+    context.active = true;
+
     GroupSelection();
     if(hover.IsEmpty() && gs.n == 0 && gs.constraints == 0) {
         // No reason to display a context menu.
@@ -591,7 +608,7 @@ void GraphicsWindow::MouseLeftDown(double mx, double my) {
 
             pending.operation = DRAGGING_NEW_LINE_POINT;
             pending.point = hr.entity(2);
-            pending.description = "click to place next point of line";
+            pending.description = "click next point of line, or press Esc";
             SK.GetEntity(pending.point)->PointForceTo(v);
             break;
 
@@ -675,7 +692,7 @@ void GraphicsWindow::MouseLeftDown(double mx, double my) {
 
             pending.operation = DRAGGING_NEW_CUBIC_POINT;
             pending.point = hr.entity(4);
-            pending.description = "click to place next point of cubic";
+            pending.description = "click next point of cubic, or press Esc";
             break;
 
         case MNU_WORKPLANE:
@@ -734,11 +751,60 @@ void GraphicsWindow::MouseLeftDown(double mx, double my) {
             break;
 
         case DRAGGING_NEW_ARC_POINT:
-        case DRAGGING_NEW_CUBIC_POINT:
             ConstrainPointByHovered(pending.point);
             ClearPending();
             break;
 
+        case DRAGGING_NEW_CUBIC_POINT: {
+            hRequest hr = pending.point.request();
+            Request *r = SK.GetRequest(hr);
+
+            if(hover.entity.v == hr.entity(1).v && r->extraPoints >= 2) {
+                // They want the endpoints coincident, which means a periodic
+                // spline instead.
+                r->type = Request::CUBIC_PERIODIC;
+                // Remove the off-curve control points, which are no longer
+                // needed here; so move [2,ep+1] down, skipping first pt.
+                int i;
+                for(i = 2; i <= r->extraPoints+1; i++) {
+                    SK.GetEntity(hr.entity((i-1)+1))->PointForceTo(
+                        SK.GetEntity(hr.entity(i+1))->PointGetNum());
+                }
+                // and move ep+3 down by two, skipping both
+                SK.GetEntity(hr.entity((r->extraPoints+1)+1))->PointForceTo(
+                  SK.GetEntity(hr.entity((r->extraPoints+3)+1))->PointGetNum());
+                r->extraPoints -= 2;
+                // And we're done.
+                SS.MarkGroupDirty(r->group);
+                SS.later.generateAll = true;
+                ClearPending();
+                break;
+            }
+
+            if(ConstrainPointByHovered(pending.point)) {
+                ClearPending();
+                break;
+            }
+
+            Entity e;
+            if(r->extraPoints >= arraylen(e.point) - 4) {
+                ClearPending();
+                break;
+            }
+
+            (SK.GetRequest(hr)->extraPoints)++;
+            SS.GenerateAll(-1, -1);
+
+            int ep = r->extraPoints;
+            Vector last = SK.GetEntity(hr.entity(3+ep))->PointGetNum();
+
+            SK.GetEntity(hr.entity(2+ep))->PointForceTo(last);
+            SK.GetEntity(hr.entity(3+ep))->PointForceTo(v);
+            SK.GetEntity(hr.entity(4+ep))->PointForceTo(v);
+            pending.point = hr.entity(4+ep);
+            break;
+        }
+
         case DRAGGING_NEW_LINE_POINT: {
             if(ConstrainPointByHovered(pending.point)) {
                 ClearPending();
@@ -755,7 +821,7 @@ void GraphicsWindow::MouseLeftDown(double mx, double my) {
             // And drag an endpoint of the new line segment
             pending.operation = DRAGGING_NEW_LINE_POINT;
             pending.point = hr.entity(2);
-            pending.description = "click to place next point of next line";
+            pending.description = "click next point of line, or press Esc";
 
             break;
         }
diff --git a/request.cpp b/request.cpp
index f45abe5..895c8ed 100644
--- a/request.cpp
+++ b/request.cpp
@@ -50,7 +50,12 @@ void Request::Generate(IdList<Entity,hEntity> *entity,
 
         case Request::CUBIC:
             et = Entity::CUBIC;
-            points = 4; 
+            points = 4 + extraPoints; 
+            break;
+
+        case Request::CUBIC_PERIODIC:
+            et = Entity::CUBIC_PERIODIC;
+            points = 3 + extraPoints;
             break;
 
         case Request::TTF_TEXT:
@@ -64,6 +69,7 @@ void Request::Generate(IdList<Entity,hEntity> *entity,
 
     // Generate the entity that's specific to this request.
     e.type = et;
+    e.extraPoints = extraPoints;
     e.group = group;
     e.style = style;
     e.workplane = workplane;
@@ -157,6 +163,7 @@ char *Request::DescriptionString(void) {
             case DATUM_POINT:       s = "datum-point"; break;
             case LINE_SEGMENT:      s = "line-segment"; break;
             case CUBIC:             s = "cubic-bezier"; break;
+            case CUBIC_PERIODIC:    s = "periodic-cubic"; break;
             case CIRCLE:            s = "circle"; break;
             case ARC_OF_CIRCLE:     s = "arc-of-circle"; break;
             case TTF_TEXT:          s = "ttf-text"; break;
diff --git a/sketch.h b/sketch.h
index e608717..7f3476c 100644
--- a/sketch.h
+++ b/sketch.h
@@ -248,11 +248,13 @@ public:
     static const int DATUM_POINT            = 101;
     static const int LINE_SEGMENT           = 200;
     static const int CUBIC                  = 300;
+    static const int CUBIC_PERIODIC         = 301;
     static const int CIRCLE                 = 400;
     static const int ARC_OF_CIRCLE          = 500;
     static const int TTF_TEXT               = 600;
 
     int         type;
+    int         extraPoints;
 
     hEntity     workplane; // or Entity::FREE_IN_3D
     hGroup      group;
@@ -302,6 +304,7 @@ public:
     static const int WORKPLANE              = 10000;
     static const int LINE_SEGMENT           = 11000;
     static const int CUBIC                  = 12000;
+    static const int CUBIC_PERIODIC         = 12001;
     static const int CIRCLE                 = 13000;
     static const int ARC_OF_CIRCLE          = 14000;
     static const int TTF_TEXT               = 15000;
@@ -313,7 +316,8 @@ public:
 
     // When it comes time to draw an entity, we look here to get the
     // defining variables.
-    hEntity     point[4];
+    hEntity     point[12];
+    int         extraPoints;
     hEntity     normal;
     hEntity     distance;
     // The only types that have their own params are points, normals,
@@ -386,6 +390,11 @@ public:
     ExprVector NormalExprsV(void);
     ExprVector NormalExprsN(void);
 
+    Vector CubicGetStartNum(void);
+    Vector CubicGetFinishNum(void);
+    ExprVector CubicGetStartTangentExprs(void);
+    ExprVector CubicGetFinishTangentExprs(void);
+
     void AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index);
     void GenerateEquations(IdList<Equation,hEquation> *l);
 };
@@ -423,6 +432,7 @@ public:
     bool IsVisible(void);
     bool PointIsFromReferences(void);
 
+    void ComputeInterpolatingSpline(SBezierList *sbl, bool periodic);
     void GenerateBezierCurves(SBezierList *sbl);
     void GenerateEdges(SEdgeList *el, bool includingConstruction=false);
 
diff --git a/solvespace.h b/solvespace.h
index 9d615fd..9ea4a86 100644
--- a/solvespace.h
+++ b/solvespace.h
@@ -33,6 +33,9 @@ inline double WRAP_SYMMETRIC(double v, double n) {
     return v;
 }
 
+// Why is this faster than the library function?
+inline double ffabs(double v) { return (v > 0) ? v : (-v); }
+
 #define SWAP(T, a, b) do { T temp = (a); (a) = (b); (b) = temp; } while(0)
 #define ZERO(v) memset((v), 0, sizeof(*(v)))
 #define CO(v) (v).x, (v).y, (v).z
@@ -234,8 +237,8 @@ bool StringEndsIn(char *str, char *ending);
 
 class System {
 public:
-    static const int MAX_UNKNOWNS = 1000;
-    static const int MAX_DRAGGED  = 4;
+    static const int MAX_UNKNOWNS      = 1024;
+    static const int MAX_DRAGGED       = 4;
 
     EntityList                      entity;
     ParamList                       param;
diff --git a/system.cpp b/system.cpp
index 4711002..38a7b1b 100644
--- a/system.cpp
+++ b/system.cpp
@@ -171,15 +171,15 @@ bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS],
         // greater. First, find a pivot (between rows i and N-1).
         max = 0;
         for(ip = i; ip < n; ip++) {
-            if(fabs(A[ip][i]) > max) {
+            if(ffabs(A[ip][i]) > max) {
                 imax = ip;
-                max = fabs(A[ip][i]);
+                max = ffabs(A[ip][i]);
             }
         }
         // Don't give up on a singular matrix unless it's really bad; the
         // assumption code is responsible for identifying that condition,
         // so we're not responsible for reporting that error.
-        if(fabs(max) < 1e-20) return false;
+        if(ffabs(max) < 1e-20) return false;
 
         // Swap row imax with row i
         for(jp = 0; jp < n; jp++) {
@@ -201,7 +201,7 @@ bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS],
     // We've put the matrix in upper triangular form, so at this point we
     // can solve by back-substitution.
     for(i = n - 1; i >= 0; i--) {
-        if(fabs(A[i][i]) < 1e-20) return false;
+        if(ffabs(A[i][i]) < 1e-20) return false;
 
         temp = B[i];
         for(j = n - 1; j > i; j--) {
@@ -294,7 +294,7 @@ bool System::NewtonSolve(int tag) {
             if(isnan(mat.B.num[i])) {
                 return false;
             }
-            if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE) {
+            if(ffabs(mat.B.num[i]) > CONVERGE_TOLERANCE) {
                 converged = false;
                 break;
             }
@@ -483,7 +483,7 @@ int System::Solve(Group *g, int *dof, List<hConstraint> *bad,
 didnt_converge:
     SK.constraint.ClearTags();
     for(i = 0; i < eq.n; i++) {
-        if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE || isnan(mat.B.num[i])) {
+        if(ffabs(mat.B.num[i]) > CONVERGE_TOLERANCE || isnan(mat.B.num[i])) {
             // This constraint is unsatisfied.
             if(!mat.eq[i].isFromConstraint()) continue;
 
diff --git a/textwin.cpp b/textwin.cpp
index 1bc630e..dd74877 100644
--- a/textwin.cpp
+++ b/textwin.cpp
@@ -335,11 +335,22 @@ void TextWindow::DescribeSelection(void) {
                     SS.MmToString((p1.Minus(p0).Magnitude())));
                 break;
             }
+            case Entity::CUBIC_PERIODIC:
             case Entity::CUBIC:
-                Printf(false, "%FtCUBIC BEZIER CURVE%E");
-                for(i = 0; i <= 3; i++) {
+                int pts;
+                if(e->type == Entity::CUBIC_PERIODIC) {
+                    Printf(false, "%FtPERIODIC C2 CUBIC SPLINE%E");
+                    pts = (3 + e->extraPoints);
+                } else if(e->extraPoints > 0) {
+                    Printf(false, "%FtINTERPOLATING C2 CUBIC SPLINE%E");
+                    pts = (4 + e->extraPoints);
+                } else {
+                    Printf(false, "%FtCUBIC BEZIER CURVE%E");
+                    pts = 4;
+                }
+                for(i = 0; i < pts; i++) {
                     p = SK.GetEntity(e->point[i])->PointGetNum();
-                    Printf((i==0), "   p%c = " PT_AS_STR, '0'+i, COSTR(p));
+                    Printf((i==0), "   p%d = " PT_AS_STR, i, COSTR(p));
                 }
                 break;
 
diff --git a/util.cpp b/util.cpp
index fa76bcf..d0a2901 100644
--- a/util.cpp
+++ b/util.cpp
@@ -125,6 +125,46 @@ void MakeMatrix(double *mat, double a11, double a12, double a13, double a14,
     mat[15] = a44;
 }
 
+//-----------------------------------------------------------------------------
+// Solve a mostly banded matrix. In a given row, there are LEFT_OF_DIAG
+// elements to the left of the diagonal element, and RIGHT_OF_DIAG elements to
+// the right (so that the total band width is LEFT_OF_DIAG + RIGHT_OF_DIAG + 1).
+// There also may be elements in the last two columns of any row. We solve
+// without pivoting.
+//-----------------------------------------------------------------------------
+void BandedMatrix::Solve(void) {
+    int i, ip, j, jp;
+    double temp;
+
+    // Reduce the matrix to upper triangular form.
+    for(i = 0; i < n; i++) {
+        for(ip = i+1; ip < n && ip <= (i + LEFT_OF_DIAG); ip++) {
+            temp = A[ip][i]/A[i][i];
+
+            for(jp = i; jp < (n - 2) && jp <= (i + RIGHT_OF_DIAG); jp++) {
+                A[ip][jp] -= temp*(A[i][jp]);
+            }
+            A[ip][n-2] -= temp*(A[i][n-2]);
+            A[ip][n-1] -= temp*(A[i][n-1]);
+
+            B[ip] -= temp*B[i];
+        }
+    }
+
+    // And back-substitute.
+    for(i = n - 1; i >= 0; i--) {
+        temp = B[i];
+
+        if(i < n-1) temp -= X[n-1]*A[i][n-1];
+        if(i < n-2) temp -= X[n-2]*A[i][n-2];
+
+        for(j = min(n - 3, i + RIGHT_OF_DIAG); j > i; j--) {
+            temp -= X[j]*A[i][j];
+        }
+        X[i] = temp / A[i][i];
+    }
+}
+
 const Quaternion Quaternion::IDENTITY = { 1, 0, 0, 0 };
 
 Quaternion Quaternion::From(double w, double vx, double vy, double vz) {
diff --git a/wishlist.txt b/wishlist.txt
index 04cca80..6081c8e 100644
--- a/wishlist.txt
+++ b/wishlist.txt
@@ -1,7 +1,6 @@
 
 multi-drag
 copy and paste
-interpolating splines
 associative entities from solid model, as a special group
 
 -----