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');
}
(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);
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);
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);
}
}
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');