+
+static void interpolate_lin4point(void) {
+ /* four points P Q R S, although P and S may be missing
+ * interpolate in QR finding M. */
+ int xq,yq,eqr, vp,vq,vm,vr,vs, k;
+ double pqtarg[D3], srtarg[D3];
+
+ if (inc != 2) fail("cannot do >2x lin4point interpolation\n");
+
+ for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */
+ printf("eqr=%d\n",eqr);
+ for (yq=0; yq<Y; yq+=inc) {
+ printf(" yq=%2d ",yq);
+ for (xq=((yq>>1)&1); xq<X; xq+=inc) {
+ vq= yq << YSHIFT | xq;
+ Traverse trav; trav.v=vq; trav.e=eqr;
+
+ traverse_next(&trav);
+ vm= trav.v;
+ if (vm<0) continue;
+
+ traverse_next(&trav);
+ vr= trav.v;
+ assert(vr>=0);
+
+ traverse_next(&trav);
+ traverse_next(&trav);
+ vs= trav.v;
+
+ trav.v= vq; trav.e= EDGE_OPPOSITE(eqr);
+ traverse_next(&trav);
+ traverse_next(&trav);
+ vp= trav.v;
+
+ printf(" 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
+ vp&0xff,vq,vm,vr,vs&0xff);
+
+ if (vp>=0)
+ K pqtarg[k]= all.a[vq][k]*2 - all.a[vp][k];
+ else
+ K pqtarg[k]= all.a[vr][k];
+
+ if (vs>=0)
+ K srtarg[k]= all.a[vr][k]*2 - all.a[vs][k];
+ else
+ K srtarg[k]= all.a[vq][k];
+
+ K {
+ const double alpha= 6;
+ all.a[vm][k]= 0.5 * (all.a[vq][k] + all.a[vr][k]);
+ pqtarg[k]= 0.5 * (pqtarg[k] + all.a[vm][k]);
+ srtarg[k]= 0.5 * (srtarg[k] + all.a[vm][k]);
+ all.a[vm][k]= (pqtarg[k] + srtarg[k] + alpha * all.a[vm][k])
+ / (2 + alpha);
+ note_computed(vm,k);
+ }
+ }
+ putchar('\n');
+ }
+ }
+}
+
+static void interpolate_gsl_line
+ (int startvertex, int direction /* edge number */, int nmax,
+ const gsl_interp_type *aperiodic, const gsl_interp_type *periodic) {