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

typedef struct {GLfloat a[3];} vec;
static vec vsub(vec a, vec b){vec d; int j=3; while(j--) d.a[j] = a.a[j] - b.a[j]; return d;}
static vec sp(double a, vec b){vec d; int j=3; while(j--) d.a[j] = a*b.a[j]; return d;} 
static vec cp(vec a, vec b){return (vec){
   {a.a[1]*b.a[2]-b.a[1]*a.a[2], a.a[2]*b.a[0]-b.a[2]*a.a[0], a.a[0]*b.a[1]-b.a[0]*a.a[1]}};}
static vec z(float * z){return (vec){{z[0], z[1], z[2]}};}

typedef unsigned char uchar;
#include "SurfaceGeometry.h"
#define num_vertices 18
typedef uchar v4[14][4];
typedef uchar v2[19][2];
typedef struct{v4 * fac; int fc; v2 * edg; int ec; uchar st;} shapE;
int const ibias = 40; // added to shapes[][].st to make invisible.
int const bias = 5; // added to shapes[][].st to enable negative offsets.
shapE shapes[4][3];

GLfloat vertices [num_vertices][3] = {
{0, 2, 4},   {1, 3, 3},    {2, 2, 2},    {1, 1, 3},    {0, 2, 2},   {1, 1, 1},
{-4, -2, 0}, {-3, -1, -1}, {-2, -2, -2}, {-3, -3, -1}, {-2, -2, 0}, {-1, -1, -1},
{-1, -1, 1}, {-3, -1, 1}, {-2, 0, 2}, {-2, -2, 2}, {-3, -1, 3}, {0, 0, 2}};

uchar faces [][4] = {
{0, 1, 2, 3}, {0, 3, 9, 6}, {0, 6, 7, 1}, {4, 3, 5, -1}, {3, 2, 5, -1}, {2, 1, 4, 5},
{3, 4, 7, 10}, {9, 8, 7, 6}, {9, 10, 11, 8}, {8, 11, 7, -1}, {11, 10, 7, -1}
// {1, 7, 4}, {9, 3, 10} empty faces to finish topological complex.
};
uchar edges [19][2] = {
{0,1}, {1,2}, {2,3}, {3,0}, {3,9}, {9,6}, {6,0}, {6,7}, {7,1}, {4,3}, {3,5},
{5,4}, {2,5}, {9,8}, {8,7}, {9,10}, {10,11}, {11,8}, {11,7}};

uchar keyFaces[][4] = {{0, 1, 3, -1}, {0, 3, 9, 6}, {0, 6, 7, 1}, {3, 1, 7, 10},
{9, 8, 7, 6}, {9, 10, 11, 8}, {8, 11, 7, -1}, {11, 10, 7, -1}};
// {9, 3, 10} empty face to finish topological complex.
uchar keyEdges[16][2] = {
{0,1}, {1,3}, {0,3}, {3,9}, {9,6}, {6,0}, {6,7}, {7,1}, {7,10}, {10,3},
{9,8}, {8,7}, {9,10}, {10,11}, {11,8}, {11,7}};

uchar augFaces[16][4] = {
{0, 1, 2, 3}, {0, 3, 12, 14}, {0, 6, 7, 1}, {4, 3, 5, -1}, {3, 2, 5, -1}, {2, 1, 4, 5},
{3, 4, 7, 10}, {9, 8, 7, 6}, {9, 10, 11, 8}, {8, 11, 7, -1}, {11, 10, 7, -1},
{9, 6, 13, 12}, {13, 15, 12, -1}, {13, 16, 15, -1}, {12, 15, 16, 14}, {14, 16, 13, -1}};
// {9, 3, 10} empty face to finish topological complex.
// {0, 14, 13, 6}

uchar augEdges[27][2] = {
{0,1}, {1,2}, {2,3}, {3,0}, {3,9}, {9,6}, {6,7}, {7,1}, {4,3}, {3,5},
{5,4}, {2,5}, {9,8}, {8,7}, {9,10}, {10,11}, {11,8}, {11,7},
{6,13}, {14,0}, {13,12}, {12,14}, {14,16}, {16,15}, {15,12}, {13,15}, {13,16}
};

static void BuildPuzzle(GLuint * polyList, GLuint * lineList){
static GLfloat woodColors [4][3] =
{{1, 0, .4}, {.3, .6, .3}, {.6, .1, .2}, {.3, .1, .6}};
	*polyList = glGenLists (1);
	glNewList(*polyList, GL_COMPILE);
	{int j = 3; while(j--) {
	{int k = 4; while(k--) {shapE sh = shapes[k][j];
	v4 * fx = sh.fac;
	void dn(int f){vec n = cp(
	         vsub(z(vertices[(*fx)[f][1]]), z(vertices[(*fx)[f][0]])),
		     vsub(z(vertices[(*fx)[f][1]]), z(vertices[(*fx)[f][2]])));
			glNormal3fv(n.a);}
	glColor3fv(woodColors[k]);
    if(sh.st < ibias/2){int K=sh.fc;
	  	while(K--){int j=4; glBegin(GL_POLYGON); dn(K);
		  while(j--) if((*fx)[K][j]<99)
		 {vec v = vsub(z(vertices[(*fx)[K][j]]), sp(sh.st-bias, z(vertices[5])));
			   glVertex3fv(v.a);}
		  glEnd ();}
		glColor3f (0, 0, 0); glBegin(GL_LINES);
		{int n = sh.ec; v2 * ex = sh.edg; while(n--) 
		{glVertex3fv(vertices[(*ex)[n][0]]); glVertex3fv(vertices[(*ex)[n][1]]);}}
		glEnd ();}

		{int k=num_vertices; while(k--) {GLfloat t = vertices[k][0];
		  vertices[k][0] = vertices[k][1]; vertices[k][1] = -t;}}}}

		{int k=num_vertices; while(k--) {GLfloat t = vertices[k][0]; vertices[k][0] = vertices[k][1];
		  vertices[k][1] = vertices[k][2]; vertices[k][2] = t;}}
		  }}
	glEndList ();
}

int pc;
int pcl;
typedef struct{uchar op; uchar opand1; uchar opand2;} inst;
inst prog[] = {
{1, 0, 2}, {1, 3, 2}, {1, 2, 0}, {1, 1, 2}, {1, 2, 2},
{1, 1, 0}, {3, 0, 2}, {3, 0, 2}, {2, 2, 0}, {1, 0, 1},
{3, 2, 0}, {2, 0, 2}, {2, 0, 2}, {3, 1, 2}, {3, 1, 2}, 
{1, 3, 0}, {3, 0, 1}, {3, 0, 1},
{2, 3, 0}, {2, 3, 0}, {1, 1, 1}, {3, 1, 0}, {3, 3, 0},
{3, 1, 1}, {2, 0, 1}, {3, 1, 0}, {2, 0, 1}, {2, 1, 1},
{2, 1, 1}, {2, 3, 0}, {1, 2, 1}, {1, 3, 1}, {3, 2, 1},
{2, 3, 1}, {2, 1, 0}, {2, 1, 0}, {3, 3, 1}, {3, 3, 1},
{3, 2, 1}, {3, 3, 0}, {3, 3, 0}, {2, 1, 2}, {2, 1, 2},
{1, 0, 0}, {3, 0, 0}, {3, 0, 0}, {3, 0, 0}};

void disassem(){
  if(!pc) return;
  inst i = prog[--pc];
  switch(i.op){ 
    case 1:  // Remove a piece
	  shapes[i.opand1][i.opand2].st += ibias; break;
	case 2: // displace -1 along body diagonal
	  --shapes[i.opand1][i.opand2].st; break;
	case 3: // displace +1 along body diagonal
	  ++shapes[i.opand1][i.opand2].st; break;
	}
return;}
void assemble(){
  if(pc == pcl) return;
  inst i = prog[pc++];
  switch(i.op){
    case 1: // Insert a piece
      shapes[i.opand1][i.opand2].st -= ibias; break;
	case 2: // displace +1 along body diagonal
	  ++shapes[i.opand1][i.opand2].st; break;
	case 3: // displace -1 along body diagonal
	  --shapes[i.opand1][i.opand2].st; break;
  }
return;}

void recompute(GLuint * polyList, GLuint * lineList){
    if (*polyList) glDeleteLists (*polyList, 1);
	if (*lineList)
	 glDeleteLists (*lineList, 1);
	*polyList = *lineList = 0;
	
	 BuildPuzzle (polyList, lineList);
}
	
void BuildGeometry (GLuint * polyList, GLuint * lineList)
{
	{int j=4*3; while(j--) {shapes[0][j] = (shapE){
	   faces, sizeof(faces)/sizeof(faces[0]),
	   edges, sizeof(edges)/sizeof(edges[0]), bias};
	   }};
	shapes[0][0] = (shapE){&keyFaces, sizeof(keyFaces)/sizeof(keyFaces[0]),
	          &keyEdges, sizeof(keyEdges)/sizeof(keyEdges[0]), bias};
    shapes[2][2] = (shapE){&augFaces, sizeof(augFaces)/sizeof(augFaces[0]),
	          &augEdges, sizeof(augEdges)/sizeof(augEdges[0]), bias};
	pc = pcl = sizeof(prog)/sizeof(prog[0]);
	    // delete existing list
recompute(polyList, lineList);}