String stp;        // used to print why we stop
float Vt=1.0, Bs=1.0;
void sfairInit() {
   laplace(L,P); laplace(B,L); mul(B,-2); mul(Q,0); mul(V,0); 
    float s=0; for (int i=0; i<n; i++) if (!F[i]) s+=dot(B[i],B[i]); s/=2;  Bs=1./s;         // early break (residues are zero)
   }

boolean sfair() {
    addScaledIfNot(P,Vt,V,F);                 // P=P+V/t   move
    addScaledIfNot(V,Bs,B,F);                 // V=V+B/s   accelerate
    laplace(L,V); laplace(Q,L); mul(Q,2.);    // Q=-3*L2(V)
    untuckIf(Q,V,F);                          // if (F) add to Q of neighbors: 2V, -0.5V, -0.5V, 2V
    float t=0; for (int i=0; i<n; i++) if (!F[i]) {t+=dot(Q[i],V[i]);}; t/=2;   Vt=1./t;  // t=sum QV
       if (Bs*Vt > 1000./6) {  stp="t"; return(false);   };  
    addScaledIfNot(B,-Vt,Q,F);               // B=B-Q/t   update acceleration to even out speeds
    float s=0; for (int i=0; i<n; i++) if (!F[i]) s+=dot(B[i],B[i]); s/=2;  Bs=1./s; 
       if (s<0.0001) { stp="s"; return(false);   };        // early break (residues are zero)
   stp=""; return(true);
    }
void drawV() {for (int i=0; i<n; i++) if(!F[i]) {vec W = V[i].make(); W.mul(Vt); W.show(P[i]);};};
void drawB() {for (int i=0; i<n; i++) if(!F[i]) {pt M = P[i].make(); vec W = V[i].make(); W.mul(Vt); M.addVec(W); W = B[i].make(); W.mul(Bs*Vt);  W.show(M);};};

void laplace(vec L[], pt P[]) {for (int i=0; i<n; i++) {L[i].setTo(P[i].vecToMid(P[pn(i)],P[nn(i)]));}; }
void laplace(vec LL[], vec L[]) {for (int i=0; i<n; i++) {LL[i].setTo(L[i].vecToMid(L[pn(i)],L[nn(i)]));}; }
void mul(vec V[], float m) {for (int i=0; i<n; i++) {V[i].mul(m);}; }
void addScaledIfNot(vec A[], float m, vec B[],boolean F[]) {for (int i=0; i<n; i++) if (!F[i]) A[i].addScaled(m,B[i]); }
void addScaledIfNot(pt A[], float m, vec B[],boolean F[]) {for (int i=0; i<n; i++) if (!F[i]) A[i].addScaledVec(m,B[i]); }
void untuckIf(vec Q[], vec V[], boolean F[]) { for (int i=0; i<n; i++) if (F[i]) { 
  Q[pn(i)].addScaled(2,V[i]); Q[nn(i)].addScaled(2,V[i]);  Q[pn(pn(i))].addScaled(-0.5,V[i]); Q[nn(nn(i))].addScaled(-0.5,V[i]); }; }

 // P = positions
 // V = velocities
 // B = accelerations

 // INIT: B = 3*L2(P),
 //            where L2(P) = (-P.p.p+4*P.p.x+4*P.n-P.n.n)/6-P = cubic fit - P = correction vector from cubic fit to P
 //            Let .n and .p denote the next and previous
 // we compute L2(P) as (-2/3)L(L(P)), where L(P) = (P.p-2P+P.n)/2

 // REPEAT:
 // s=sum B^2
 // V=V+B/s   accelerate
 // Q=-3*L2(V)
 // if (F) add to Q of neighbors: 2V, -0.5V, -0.5V, 2V
 // t=sum QV
 // B=B-Q/t   update acceleration to even out speeds
 // P=P+V/t   move