#include "mgraph.h"
-/* filled in by characterise_input: */
-static int shift, oldxbits, oldybits, oldx, oldy, oldsz, inc;
+#define OXBITS (XBITS-1)
+#define OX (1<<OXBITS)
+#define OY (Y/2 + 1)
+#define ON (OX*OY)
+#define INC 2
/* filled in by read_input: */
static struct Vertices all;
+static int computed_count[N];
static void characterise_input(void) {
struct stat stab;
if (!stab.st_size || stab.st_size % (sizeof(double)*D3) ||
stab.st_size > INT_MAX)
- fail("input file is not reasonable whole number of doubles\n");
-
- oldsz= stab.st_size / (sizeof(double)*D3);
- for (shift=1;
- shift < XBITS+1 && shift < YBITS+1;
- shift++) {
- oldxbits= XBITS-1;
- oldybits= YBITS-1;
- oldx= 1<<oldxbits;
- oldy= (1<<oldybits)-1;
- fprintf(stderr,"sizeof(double)=%d XYBITS=%d,%d, XY=%d*%d=%d"
- " oldsz=%d shift=%d oldxybits=%d,%d oldxy=%d*%d=%d\n",
- (int)sizeof(double), XBITS,YBITS, X,Y,N,
- oldsz,shift,oldxbits,oldybits,oldx,oldy,oldx*oldy);
- if (oldx*oldy == oldsz) goto found;
- }
- fail("input file size cannot be interpolated to target file size\n");
+ fail("input file is not reasonable whole number of vertices\n");
+
+ fprintf(stderr,"XBITS=%d XY=%d*%d=%d\n", XBITS, X,Y,N);
+ fprintf(stderr,"OXBITS=%d OXY=%d*%d=%d\n", OXBITS, OX,OY,ON);
+ fprintf(stderr,"sizeof(double)=%d st_size=%lu\n",
+ (int)sizeof(double), (unsigned long)stab.st_size);
- found:
- inc= 1<<shift;
- fprintf(stderr,"inc=%d\n",inc);
+ if (stab.st_size != (ON*sizeof(double)*D3))
+ fail("input file wrong size\n");
}
static void read_input(void) {
int x,y, ox,oy, v,ov;
- for (oy=y=0; y<Y; oy++, y+=inc) {
+ for (oy=y=0; y<Y; oy++, y+=INC) {
fprintf(stderr, "y=%2d>%2d", oy,y);
- for (ox=x=0; x<X; ox++, x+=inc) {
+ for (ox=x=0; x<X; ox++, x+=INC) {
errno= 0;
- ov= (oy << oldxbits) | ox;
+ ov= (oy << OXBITS) | ox;
v= (y << YSHIFT) | x;
fprintf(stderr, " 0%02o->0x%02x", ov, v);
if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
diee("\nread input");
+ computed_count[v]= -1;
}
fputc('\n',stderr);
}
*
*/
-#define NEXTV (v= EDGE_END2(v,direction))
-
+typedef struct {
+ int v, e;
+} Traverse;
+
+static void traverse_next(Traverse *t) {
+ int v2;
+ v2= EDGE_END2(t->v, t->e);
+ if (v2>=0) {
+ int e2= edge_reverse(t->v, t->e);
+ assert(EDGE_END2(v2,e2) == t->v);
+ t->e= (e2+3) % V6;
+ }
+ t->v= v2;
+}
+
static void interpolate_line(int startvertex,
int direction /* edge number */,
- int norig) {
- double xa[norig+1], ya[D3][norig+1], n;
+ int nmax) {
+ double xa[nmax], ya[D3][nmax];
const gsl_interp_type *it;
gsl_interp *interp;
gsl_interp_accel *accel;
- int i, v, k;
+ int i, k, nold, nnew;
+ Traverse traverse;
+
+#define STARTV (traverse.v=startvertex, traverse.e=direction)
+#define NEXTV (traverse_next(&traverse), fprintf(stderr,"-%02x",traverse.v))
+
+ fprintf(stderr,"interpolate_line 0x%02x,%d nmax=%d ",
+ startvertex,direction,nmax);
- for (i=0, v=startvertex;
+ for (i=0, STARTV;
;
- i++, NEXTV, assert(v>=0), NEXTV) {
- if (v<0)
- break;
- assert(i <= norig);
+ ) {
+ assert(traverse.v>=0);
+ assert(i < nmax);
xa[i]= i;
- K ya[k][i]= all.a[v][k];
- if (v==startvertex)
- break;
+ K ya[k][i]= all.a[traverse.v][k];
+ fputc('*',stderr);
+ if (i && traverse.v==startvertex) { fputc('@',stderr); break; }
+ i++;
+ NEXTV;
+ if (traverse.v<0) break;
+ NEXTV;
}
- n= i;
- it= v>=0 ? gsl_interp_cspline_periodic : gsl_interp_cspline;
+ if (traverse.v>=0) {
+ it= gsl_interp_cspline_periodic;
+ nold= i+1;
+ nnew= i;
+ } else {
+ it= gsl_interp_cspline;
+ nold= i;
+ nnew= i-1;
+ }
+
+ fprintf(stderr," n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
K {
- interp= gsl_interp_alloc(it,n); if (!interp) diee("gsl_interp_alloc");
+ fprintf(stderr,"\n k=%d ", k);
+
+ interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc");
accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
- GA( gsl_interp_init(interp,xa,ya[k],n) );
+ GA( gsl_interp_init(interp,xa,ya[k],nold) );
- for (i=0, v=startvertex;
- i<n-1;
+ for (i=0, STARTV;
+ i<nnew;
i++, NEXTV) {
- assert(v>=0); NEXTV; assert(v>=0);
- GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel, &all.a[v][k]) );
+ assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
+ fputc('#',stderr);
+ assert(computed_count[traverse.v] == k);
+ GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
+ &all.a[traverse.v][k]) );
+ computed_count[traverse.v]++;
}
gsl_interp_accel_free(accel);
gsl_interp_free(interp);
}
+ fprintf(stderr," done\n");
}
static void interpolate(void) {
int x,y;
- if (shift!=1) fail("only interpolation by factor of 2 supported\n");
-
- for (y=0; y<Y; y+=inc) {
- interpolate_line(y<<YSHIFT, 0, oldx);
+ for (y=0; y<(Y+1)/2; y+=INC) {
+ interpolate_line(y<<YSHIFT, 0, OX*2+1);
}
- for (x=0; x<X; x+=inc) {
- interpolate_line(x, 5, oldy);
- interpolate_line(x, 4, oldy);
+ for (x=0; x<X; x+=INC) {
+ interpolate_line( x, 5, OY);
+ interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
}
}
static void write_output(void) {
+ int x, y, v, c, bad=0;
+ for (y=0; y<Y; y++) {
+ fprintf(stderr,"checking y=%d ", y);
+ for (x=0; x<X; x++) {
+ v= y<<YSHIFT | x;
+ fprintf(stderr," %02x",v);
+ c= computed_count[v];
+ if (c==D3) fputc('#',stderr);
+ else if (c==-1) fputc('*',stderr);
+ else { fputc('!',stderr); bad++; }
+ }
+ fputc('\n',stderr);
+ }
+ assert(!bad);
if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
+ fprintf(stderr,"interpolated\n");
}
int main(int argc, const char **argv) {
~ian / moebius2.git / blobdiff