another rune

[moebius2.git] / interpolate.c

diff --git a/interpolate.c b/interpolate.c
index 9225b8f3ef5e6e241d83b5463cdb366ec0373a25..a2d1cafd48de65ad50b943f92afd6495a11de0c7 100644 (file)
--- a/interpolate.c
+++ b/interpolate.c
@@ -56,18 +56,18 @@ static void read_input(void) {
   int x,y, ox,oy, v,ov;
   
   for (oy=0; oy<oY; oy++) {
-    y= oy*2;
+    y= oy*inc;
     printf("y=%2d>%2d", oy,y);
     for (ox=0; ox<oX; ox++) {
-      x= ox*2 + (oy&1);
-      if (x>=X) { y= (Y-1)-y; x-=X; }
+      x= (ox - (oy>>1))*inc + (y>>1) + X*2;
+      while (x>=X) { y= (Y-1)-y; x-=X; }
       errno= 0;
       ov= (oy << oXBITS) | ox;
       v= (y << YSHIFT) | x;
       printf(" 0%02o->0x%02x", ov, v);
       if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
        diee("\nread input");
-      computed_count[v]= -1;
+      computed_count[v]= D3*2;
     }
     putchar('\n');
   }
@@ -173,43 +173,62 @@ static void interpolate_gsl_line
   (int startvertex, int direction /* edge number */, int nmax,
    const gsl_interp_type *aperiodic, const gsl_interp_type *periodic) {
   double xa[nmax], ya[D3][nmax];
+  int vnew[nmax], xnew[nmax];
   const gsl_interp_type *it;
   gsl_interp *interp;
   gsl_interp_accel *accel;
-  int i, k, nold, nnew;
+  int i, k, nold, nnew, x, skipped, realstart;
   Traverse traverse;
 
-#define STARTV (traverse.v=startvertex, traverse.e=direction)
-#define NEXTV (traverse_next(&traverse), printf("-%02x",traverse.v))
-
   printf("interpolate_line 0x%02x,%d nmax=%d ",
          startvertex,direction,nmax);
 
-  for (i=0, STARTV;
-       ;
-       ) {
-    assert(traverse.v>=0);
-    assert(i < nmax);
-    xa[i]= i;
-    K ya[k][i]= all.a[traverse.v][k];
-    putchar('*');
-    if (i && traverse.v==startvertex) { putchar('@'); break; }
-    i++;
-    NEXTV;
-    if (traverse.v<0) break;
-    NEXTV;
+  traverse.v=startvertex, traverse.e=direction;
+  skipped= 0;
+  while (!computed_count[traverse.v]) {
+    printf("-%02x~",traverse.v);
+    traverse_next(&traverse);
+    skipped++;
+    assert(skipped < nmax);
   }
-  if (traverse.v>=0) {
-    it= periodic;
-    nold= i+1;
-    nnew= i;
-  } else {
-    it= aperiodic;
-    nold= i;
-    nnew= i-1;
+  realstart= traverse.v;
+
+  nold=0, nnew=0;
+  for (x=0;
+       ;
+       x++, traverse_next(&traverse)) {
+    printf("-%02x",traverse.v);
+    if (traverse.v<0) {
+      putchar('$');
+      it= aperiodic;
+      assert(!skipped);
+      break;
+    }
+    if (computed_count[traverse.v]) {
+      assert(nold < nmax);
+      xa[nold]= x;
+      K ya[k][nold]= all.a[traverse.v][k];
+      nold++;
+      putchar('*');
+    } else {
+      assert(nnew < nmax);
+      xnew[nnew]= x;
+      vnew[nnew]= traverse.v;
+      nnew++;
+      putchar('#');
+    }
+    if (x && traverse.v==realstart) {
+      putchar('@');
+      it= periodic;
+      break;
+    }
   }
 
-  printf("  n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
+  printf("  n=%d->%d", nold,nnew);
+  if (!nnew) { printf(" nothing\n"); return; }
+
+  assert(nold);
+  printf(" x=%g..%g->%d..%d", xa[0],xa[nold-1], xnew[0],xnew[nnew-1]);
 
   K {
     printf("\n  k=%d ", k);
@@ -218,14 +237,12 @@ static void interpolate_gsl_line
     accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
     GA( gsl_interp_init(interp,xa,ya[k],nold) );
 
-    for (i=0, STARTV;
-        i<nnew;
-        i++, NEXTV) {
-      assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
-      putchar('#');
-      GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
-                           &all.a[traverse.v][k]) );
-      note_computed(traverse.v, k);
+    for (i=0; i<nnew; i++) {
+      int v= vnew[i];
+      x= xnew[i];
+      printf(" %02x(%d)",v,x);
+      GA( gsl_interp_eval_e(interp,xa,ya[k], x, accel, &all.a[v][k]) );
+      note_computed(v,k);
     }
     
     gsl_interp_accel_free(accel);
@@ -238,13 +255,13 @@ static void interpolate_gsl(const gsl_interp_type *aperiodic,
                            const gsl_interp_type *periodic) {
   int x,y;
   
-  for (y=0; y<(Y+1)/2; y+=inc) {
-    interpolate_gsl_line(y<<YSHIFT | ((y>>1)&1), 0, oX*2+1,
-                        aperiodic, periodic);
-  }
   for (x=0; x<X; x+=inc) {
-    interpolate_gsl_line(                x, 5, oY, aperiodic, periodic);
-    interpolate_gsl_line((Y-1)<<YSHIFT | x, 1, oY, aperiodic, periodic);
+    interpolate_gsl_line(                x, 5, Y, aperiodic, periodic);
+    interpolate_gsl_line((Y-1)<<YSHIFT | x, 1, Y, aperiodic, periodic);
+  }
+  for (y=0; y<(Y+1)/2; y++) {
+    interpolate_gsl_line(y<<YSHIFT | ((y>>1)&1), 0, X*2+1,
+                        aperiodic, periodic);
   }
 }
 
@@ -260,7 +277,7 @@ static void write_output(const char *fn) {
       printf(" %02x",v);
       c= computed_count[v];
       if (c==D3) putchar('#');
-      else if (c==-1) putchar('*');
+      else if (c==D3*2) putchar('*');
       else { printf("!%d",c); bad++; }
     }
     putchar('\n');