/*
 *   icx =  -2   -1    0   +1   +2   +3
 *          \/   \/   \/   \/   \/   \/
 *  ...___________________________________...   triangle    point
 *     /\        /\        /\        /\
 *    /  \      /  \      /  \      /  \      .
 *   /    \    /    \    /    \    /    \    .   icy = 2
 *  /      \  /      \  /      \  /      \  /
 * /________\/________\/________\/________\/                icy = 2
 * \        /\        /\        /\        /\
 *  \      /  \      /  \      /  \      /  \
 *   \    /    \    /    \    /    \    /    .   icy = 1
 *    \  /      \  /      \  /      \  /      .
 *     \/________\/________\/________\/____...              icy = 1
 *     /\        /\        /\        /\
 *    /  \      /  \      /  \      /  \      .
 *   /    \    /    \    /    \    /    \    .   icy = 0
 *  /      \  /      \  /      \  /      \  /
 * /________\/________\/________\/________\/                icy = 0
 * \        /\        /\        /\        /\
 *  \      /  \      /  \      /  \      /  \
 *   \    /    \    /    \    /    \    /    .   icy = -1
 *    \  /      \  /      \  /      \  /      .
 *     \/________\/________\/________\/____...              icy = -1
 *     /\        /\        /\        /\
 *    /  \      /  \      /  \      /  \      .
 *   /    \    /    \    /    \    /    \    .   icy = -2
 *  /      \  /      \  /      \  /      \  /
 * /________\/________\/________\/________\/
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>

#define SQRT_3_4 .86602540378443864676372317 /* sqrt(3/4) */

typedef struct pt PT;
typedef struct ln LN;
typedef struct tri TRI;
typedef struct list LIST;
typedef struct shape SHAPE;
typedef struct ll LL;

#define LPL 10

struct pt {
  double x;
  double y;
  int ix;
  int iy;
  } ;

struct ln {
  int pts[2];
  } ;

struct tri {
  int pts[3];
  int lns[3];
  int ngs[3];
  double cx;
  double cy;
  int icx;
  int icy;
  } ;

struct list {
  unsigned int nlev;
  unsigned int size;
  unsigned int maxsize;
  LL *v;
  } ;

struct ll {
  void *v[LPL];
  } ;

struct shape {
  int nt;
  int *t;
  int dup;
  } ;

static LIST pts;
static LIST lns;
static LIST trs;
static LIST shs;

static void list_init(LIST *l)
{
 l->nlev = 0;
 l->size = 0;
 l->maxsize = LPL;
 l->v = malloc(sizeof(LL));
}

static void list_walk(LIST *l, void (*fn)(unsigned int, void *))
{
 __label__ done;

 int x;

 void walk(LL *ll, unsigned int nl)
  { int i;
    for (i=0;i<LPL;i++)
     { if (nl < 1)
	{ if (x >= l->size) goto done;
	  (*fn)(x,ll->v[i]);
	  x ++;
	}
       else
	{ walk(ll->v[i],nl-1);
	}
     }
  }

 x = 0;
 walk(l->v,l->nlev);
done:;
}

static void **list_loc(LIST *l, unsigned int x, int create)
{
 LL *ll;
 void *vp;
 unsigned int nl;
 unsigned int d;
 int i;

 ll = l->v;
 d = l->maxsize / LPL;
 nl = l->nlev;
 while (nl > 0)
  { i = x / d;
    x %= d;
    vp = ll->v[i];
    if (! vp)
     { LL *new;
       if (! create) abort();
       new = malloc(sizeof(LL));
       ll->v[i] = new;
       vp = new;
     }
    ll = vp;
  }
 return(&ll->v[x]);
}

static void list_add(LIST *l, void *item)
{
 if (l->size >= l->maxsize)
  { LL *ll;
    ll = malloc(sizeof(LL));
    ll->v[0] = l->v;
    l->v = ll;
    l->maxsize *= LPL;
    l->nlev ++;
  }
 *list_loc(l,l->size++,1) = item;
}

static void *list_get(LIST *l, unsigned int x)
{
 return(*list_loc(l,x,0));
}

static unsigned int list_size(LIST *l)
{
 return(l->size);
}

static void find_tri(int icx, int icy)
{
 XXX XXX XXX;;
}

static void init(void)
{
 list_init(&pts);
 list_init(&lns);
 list_init(&trs);
 list_init(&shs);
 find_tri(0,0);
}

static void gen_shapes(void)
{
 int size;
 int prevmin;
 int prevmax;
 int curmin;
 int curmax;
 int sx;
 SHAPE s;
 int tx;
 int nx;
 int nt;
 SHAPE ns;
 int i;

 shapes[0] = 1ULL << MID_TRI;
 nshapes = 1;
 curmin = 0;
 curmax = 0;
 for (size=2;size<=MAX_SIZE;size++)
  { prevmin = curmin;
    prevmax = curmax;
    /* printf("size=%d (prev min=%d max=%d)\n",size,prevmin,prevmax); */
    curmin = MAX_SHAPES;
    curmax = -1;
    for (sx=prevmin;sx<=prevmax;sx++)
     { s = shapes[sx];
       /*
       printf("shape ");
       print_shape(s);
       printf(" from %d\n",sx);
	*/
       for (tx=N_TRIS-1;tx>=0;tx--)
	{ if (s & (1ULL << tx))
	   { for <"ngh"> (nx=3-1;nx>=0;nx--)
	      { nt = tris[tx].ngs[nx];
		if (nt < 0) continue;
		if (s & (1ULL << nt)) continue;
		ns = s | (1ULL << nt);
		/* printf("tx=%d nx=%d nt=%d ns=",tx,nx,nt); */
		/* print_shape(ns);*/
		for (i=curmin;i<=curmax;i++)
		 { if (shapes[i] == ns)
		    { /* printf(" duplicate of %d\n",i); */
		      continue <"ngh">;
		    }
		 }
		/* printf(" saving as %d\n",nshapes); */
		if (nshapes >= MAX_SHAPES) abort();
		shapes[nshapes] = ns;
		if (nshapes < curmin) curmin = nshapes;
		if (nshapes > curmax) curmax = nshapes;
		nshapes ++;
	      }
	   }
	}
     }
    /* printf("size loop done, min=%d max=%d\n",curmin,curmax); */
  }
}

int main(int, char **);
int main(int ac, char **av)
{
 if (ac != 2) usage();
 init();
 gen_shapes();
 filter_shapes();
 return(0);
}

#if 0

static int find_ln(int p1, int p2)
{
 int i;

 for (i=N_LNS-1;i>=0;i--)
  { if ( ((lns[i].pts[0] == p1) && (lns[i].pts[1] == p2)) ||
	 ((lns[i].pts[1] == p1) && (lns[i].pts[0] == p2)) ) return(i);
  }
 abort();
}

static int find_ng(int p1, int p2, int exc)
{
 int i;

 for (i=N_TRIS-1;i>=0;i--)
  { if (i == exc) continue;
    if ( ((tris[i].pts[0] == p1) && (tris[i].pts[1] == p2)) ||
	 ((tris[i].pts[1] == p1) && (tris[i].pts[0] == p2)) ||
	 ((tris[i].pts[1] == p1) && (tris[i].pts[2] == p2)) ||
	 ((tris[i].pts[2] == p1) && (tris[i].pts[1] == p2)) ||
	 ((tris[i].pts[2] == p1) && (tris[i].pts[0] == p2)) ||
	 ((tris[i].pts[0] == p1) && (tris[i].pts[2] == p2)) ) return(i);
  }
 return(-1);
}

static int hash_tri(double cx, double cy)
{
 return((11*(3+(int)rint(cy/SQRT_3_4)))+5+(int)rint(cx/.5));
}

static void gen_mesh(void)
{
 int i;
 double d;
 int errs;
 int h;

 errs = 0;
 for (i=N_LNS-1;i>=0;i--)
  { d = hypot( pts[lns[i].pts[0]].x - pts[lns[i].pts[1]].x,
	       pts[lns[i].pts[0]].y - pts[lns[i].pts[1]].y );
    if ((d < .99) || (d > 1.01))
     { printf("ln %d: %d-%d: d=%g\n",i,lns[i].pts[0],lns[i].pts[1],d);
       errs = 1;
     }
  }
 for (i=THASHSIZE-1;i>=0;i--) thash[i] = -1;
 for (i=N_TRIS-1;i>=0;i--)
  { tris[i].lns[0] = find_ln(tris[i].pts[0],tris[i].pts[1]);
    tris[i].lns[1] = find_ln(tris[i].pts[0],tris[i].pts[2]);
    tris[i].lns[2] = find_ln(tris[i].pts[1],tris[i].pts[2]);
    tris[i].ngs[0] = find_ng(tris[i].pts[1],tris[i].pts[2],i);
    tris[i].ngs[1] = find_ng(tris[i].pts[0],tris[i].pts[2],i);
    tris[i].ngs[2] = find_ng(tris[i].pts[0],tris[i].pts[1],i);
    tris[i].cx = (pts[tris[i].pts[0]].x + pts[tris[i].pts[1]].x + pts[tris[i].pts[2]].x) / 3;
    tris[i].cy = (pts[tris[i].pts[0]].y + pts[tris[i].pts[1]].y + pts[tris[i].pts[2]].y) / 3;
    h = hash_tri(tris[i].cx,tris[i].cy);
    if ((h < 0) || (h >= THASHSIZE) || (thash[h] >= 0)) abort();
    thash[h] = i;
  }
 if (errs) exit(1);
}

static void print_shape(SHAPE s)
{
 int i;
 const char *pre;

 printf("<");
 pre = "";
 for (i=0;i<N_TRIS;i++)
  { if (s & (1ULL << i))
     { printf("%s%d",pre,i);
       pre = ",";
     }
  }
 printf(">");
}

/*
 * Compute the transformation of the plane that takes fp1 to tp1, fp2
 *  to tp2, and fp3 to tp3.  Transform s with this transformation; if
 *  the result is not equal to s, call (*fn) with the new SHAPE.
 *
 * To do this, we solve for M and b in
 *
 *	[tp1.x]   [ Mxx Mxy ] [ fp1.x ]   [b.x]
 *	[     ] = [         ] [       ] + [   ]
 *	[tp1.y]   [ Myx Myy ] [ fp1.y ]   [b.y]
 *
 *	[tp2.x]   [ Mxx Mxy ] [ fp2.x ]   [b.x]
 *	[     ] = [         ] [       ] + [   ]
 *	[tp2.y]   [ Myx Myy ] [ fp2.y ]   [b.y]
 *
 *	[tp3.x]   [ Mxx Mxy ] [ fp3.x ]   [b.x]
 *	[     ] = [         ] [       ] + [   ]
 *	[tp3.y]   [ Myx Myy ] [ fp3.y ]   [b.y]
 *
 *  This gives six equations in six unknowns:
 *
 *	[1 0 0 0 0 0] [tp1.x]   [fp1.x fp1.y   0     0     1     0  ] [Mxx]
 *	[0 1 0 0 0 0] [tp1.y]   [  0     0   fp1.x fp1.y   0     1  ] [Mxy]
 *	[0 0 1 0 0 0] [tp2.x] = [fp2.x fp2.y   0     0     1     0  ] [Myx]
 *	[0 0 0 1 0 0] [tp2.y]   [  0     0   fp2.x fp2.y   0     1  ] [Myy]
 *	[0 0 0 0 1 0] [tp3.x]   [fp3.x fp3.y   0     0     1     0  ] [b.x]
 *	[0 0 0 0 0 1] [tp3.y]   [  0     0   fp3.x fp3.y   0     1  ] [b.y]
 *
 *  We solve this with traditional Gauss-Jordan elimination.  Because
 *  of the way our arguments are generated, it is a can't-happen for
 *  the matrix on the RHS above to be anywhere near ill-determined.
 *
 * Once we have this M and b, we identify the transforms of all the
 *  triangles in the shape, by transforming the centre of the source
 *  triangle according to M and b, then using hash_tri on it.  If the
 *  result is outside thash, if the thash bucket indicated is empty, or
 *  if the centre of the triangle indicated by the thash[] bucket is
 *  more than a small distance from the transformed centre
 *  (mathematically zero, but numerical errors occur), then the
 *  transformed triangle is outside the mesh; this is a can't-happen.
 *  Otherwise, this gives us a list of triangles making up the
 *  transformed piece.  It is another can't-happen for any particular
 *  triangle to occur more than once in this list.  We construct the
 *  corresponding SHAPE; if it is different from s, we pass it to
 *  (*fn)().
 */
static void maybe_call_transformed(
	void (*fn)(SHAPE), SHAPE s,
	int fp1, int fp2, int fp3,
	int tp1, int tp2, int tp3 )
{
 double eq[12][6];
 int x;
 int y;
 int my;
 double p;
 int i;
 double t;
 double Mb[6];
 double tpv[6];
 int tx;
 double fcx;
 double fcy;
 double tcx;
 double tcy;
 int ttx;
 int hb;
 SHAPE s2;

 static void dump_eq(const char *tag __attribute__((__unused__)))
  { /*
    int xi;
    int yi;
    printf("%s\n",tag);
    for (yi=0;yi<6;yi++)
     { for (xi=0;xi<12;xi++)
	{ printf(" %10.5f",eq[xi][yi]);
	}
       printf("\n");
     }
     */
  }

 for (x=12-1;x>=0;x--) for (y=6-1;y>=0;y--) eq[x][y] = 0;
 eq[0][0] = pts[fp1].x;
 eq[2][1] = pts[fp1].x;
 eq[1][0] = pts[fp1].y;
 eq[3][1] = pts[fp1].y;
 eq[0][2] = pts[fp2].x;
 eq[2][3] = pts[fp2].x;
 eq[1][2] = pts[fp2].y;
 eq[3][3] = pts[fp2].y;
 eq[0][4] = pts[fp3].x;
 eq[2][5] = pts[fp3].x;
 eq[1][4] = pts[fp3].y;
 eq[3][5] = pts[fp3].y;
 eq[4][0] = 1;
 eq[5][1] = 1;
 eq[4][2] = 1;
 eq[5][3] = 1;
 eq[4][4] = 1;
 eq[5][5] = 1;
 for (x=6-1;x>=0;x--) eq[x+6][x] = 1;
 dump_eq("initial");
 for (x=6-1;x>=0;x--)
  { my = x;
    for (y=x-1;y>=0;y--) if (fabs(eq[x][y]) > fabs(eq[x][my])) my = y;
    p = eq[x][my];
    /* printf("x=%d my=%d p=%g\n",x,my,p); */
    if (fabs(p) < .0001) abort();
    if (my != x)
     { for (i=12-1;i>=0;i--)
	{ t = eq[i][my];
	  eq[i][my] = eq[i][x];
	  eq[i][x] = t;
	}
       dump_eq("after row swap");
     }
    for (i=12-1;i>=0;i--) eq[i][x] /= p;
    dump_eq("after dividing");
    for (y=6-1;y>=0;y--)
     { if (y != x)
	{ t = eq[x][y];
	  for (i=12-1;i>=0;i--) eq[i][y] -= t * eq[i][x];
	}
     }
    dump_eq("after subtracting");
  }
 tpv[0] = pts[tp1].x;
 tpv[1] = pts[tp1].y;
 tpv[2] = pts[tp2].x;
 tpv[3] = pts[tp2].y;
 tpv[4] = pts[tp3].x;
 tpv[5] = pts[tp3].y;
 for (y=6-1;y>=0;y--)
  { t = 0;
    for (x=6-1;x>=0;x--) t += tpv[x] * eq[x+6][y];
    Mb[y] = t;
  }
 /* printf("Mb = %g %g %g %g %g %g\n",Mb[0],Mb[1],Mb[2],Mb[3],Mb[4],Mb[5]); */
 s2 = 0;
 for (tx=N_TRIS-1;tx>=0;tx--)
  { if (s & (1ULL << tx))
     { fcx = tris[tx].cx;
       fcy = tris[tx].cy;
       tcx = (Mb[0] * fcx) + (Mb[1] * fcy) + Mb[4];
       tcy = (Mb[2] * fcx) + (Mb[3] * fcy) + Mb[5];
       hb = hash_tri(tcx,tcy);
       if ((hb < 0) || (hb >= THASHSIZE) || (thash[hb] < 0)) abort();
       ttx = thash[hb];
       if ((fabs(tris[ttx].cx-tcx) > .01) || (fabs(tris[ttx].cy-tcy) > .01)) abort();
       if (s2 & (1ULL << ttx)) abort();
       s2 |= 1ULL << ttx;
     }
  }
 if (s2 != s) (*fn)(s2);
}

static void map_varieties(SHAPE s, void (*fn)(SHAPE))
{
 int t;

 for (t=N_TRIS-1;t>=0;t--)
  { if (s & (1ULL << t))
     { maybe_call_transformed( fn, s,
		tris[t].pts[0], tris[t].pts[1], tris[t].pts[2],
		tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] );
       maybe_call_transformed( fn, s,
		tris[t].pts[0], tris[t].pts[2], tris[t].pts[1],
		tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] );
       maybe_call_transformed( fn, s,
		tris[t].pts[1], tris[t].pts[2], tris[t].pts[0],
		tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] );
       maybe_call_transformed( fn, s,
		tris[t].pts[1], tris[t].pts[0], tris[t].pts[2],
		tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] );
       maybe_call_transformed( fn, s,
		tris[t].pts[2], tris[t].pts[0], tris[t].pts[1],
		tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] );
       maybe_call_transformed( fn, s,
		tris[t].pts[2], tris[t].pts[1], tris[t].pts[0],
		tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] );
     }
  }
}

static void filter_shapes(void)
{
 int i;
 int j;
 SHAPE s;

 void drop(SHAPE s2)
  { int k;
    for (k=i+1;k<nshapes;k++)
     { if (shapes[k] == s2)
	{ shapedup[k] = 1;
	  return;
	}
     }
  }

 bzero(&shapedup[0],MAX_SHAPES);
 j = 0;
 for (i=0;i<nshapes;i++)
  { if (shapedup[i]) continue;
    s = shapes[i];
    shapes[j++] = s;
    printf("Piece: ");
    print_shape(s);
    printf("\n");
    map_varieties(s,&drop);
  }
 nshapes = j;
}

static void usage(void)
{
 fprintf(stderr,"Usage: %s max-triangles\n",__progname);
 exit(1);
}

#endif