/*usr/bin/gcc -O3 -g rubik.c -lglut -lGLU -lGL -lreadline -o rubik ; exit
 *
 *      A toy program for group-based modelling of Rubik's cube
 *      and its 3D visualization using OpenGL
 *
 *      Copyright (c) 2019 Jiri Pittner <jiri@pittnerovi.com>
 *
 *      This program is free software: you can redistribute it and/or modify
 *      it under the terms of the GNU General Public License as published by
 *      the Free Software Foundation, either version 3 of the License, or
 *      (at your option) any later version.
 *
 *      This program is distributed in the hope that it will be useful,
 *      but WITHOUT ANY WARRANTY; without even the implied warranty of
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *      GNU General Public License for more details.
 *      You should have received a copy of the GNU General Public License
 *      along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
*/

//TODO (but will probably not have time to)
//group membership by the two theorems of cube group (independently of GAP)
//solving by other published algorithms (independently of GAP)
//NxNxN cube generalization - face centers of cubicles would get coordinates transformed by layer rotations (N or N-1 x generators, y generators, z generators) and permutations of cubicle faces would be computed by comparison of the resulting coordinates with the original ones 
//In this way the group can be constructed for the general case


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <sys/wait.h>
#include <unistd.h>
#include <ctype.h>
#include <signal.h>
#include <fcntl.h>
#include <time.h>
#include <signal.h>


#include <readline/readline.h>
#include <readline/history.h>


#include <GL/gl.h>
#include <GL/glut.h>
#include <GL/freeglut_ext.h>



/*
 * NOTE: I initially decided to separately number edges and vertices in a hope for more efficiency
 * but it turned out it brings nothing and just complicates the interface to GAP and to the user
 * but I am not going to rewrite it now ...
 *
 * NOTE also that we use synonymously group elements are their action on the ordered cube
 *
 * cube numbering clockwise separately edges and vertices:
 * (is arbitrary but must be consistent between generator definition and cube_numbering array)
 * Up 1-4
 * Back 5-8
 * Front 9-12
 * Down 13-16
 * Right 17-20
 * Left 21-24
 *
 * Edges: 1-5, 2-17, 3-9, 4-21, 10-20, 8-18, 6-24, 12-22, 11-13, 14-19, 7-15, 16-23
 * Vertices: 1-24-5, 2-8-17, 3-20-9, 4-12-21, 10-19-13, 7-14-18, 6-23-15, 11-16-22
 * Generators: (edge)(edge) (vertex)(vertex)(vertex)
 * U = (1 2 3 4) (9 21 5 17) (1 2 3 4) (9 21 5 17) (12 24 8 20)
 * B = (5 6 7 8) (1 24 15 18) (5 6 7 8) (2 24 15 18) (1 23 14 17)
 * F = (9 10 11 12) (13 22 3 20) (9 10 11 12) (13 22 4 20) (16 21 3 19)
 * D = (13 14 15 16) (7 23 11 19) (13 14 15 16) (7 23 11 19) (6 22 10 18)
 * R = (17 18 19 20) (10 2 8 14) (17 18 19 20) (10 3 8 14) (9 2 7 13)
 * L = (21 22 23 24) (4 12 16 6) (21 22 23 24) (12 16 6 1) (11 15 5 4)
 *
 * Cube orientation for visualization
 * Right = +x
 * Back = +y
 * Up= +z
*/

//Rubik's group is a subgroup of S24 x S24
#define NP 24

const char faces[] = "UBFDRL";

//numbering at each face as above, while coordinates
//of the boxes are in ascending xy, xz, or yz going as 00,01,02,10,11,12,20,21,22
const int cube_numbering[6][3][3] = {
{{4,4,1},{3,0,1},{3,2,2}}, //U
{{6,6,5},{7,0,5},{7,8,8}}, //B
{{11,12,12},{11,0,9},{10,10,9}}, //F
{{16,16,15},{13,0,15},{13,14,14}}, //D
{{19,20,20},{19,0,17},{18,18,17}}, //R
{{22,22,21},{23,0,21},{23,24,24}}, //L
};

const int is_edge[3][3] = {
{0,1,0},
{1,0,1},
{0,1,0}
};

const int is_vertex[3][3] = {
{1,0,1},
{0,0,0},
{1,0,1}
};



//define formal colors of original cube state 
const int edgecolors[NP+1] = {-1,0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5};
const int vertexcolors[NP+1] = {-1,0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5};

//define coordinates of face planes (x=0,y=1,z=2)
const int face_planes_and_normal[6][3] = {
{0,1,2},
{0,2,1},
{0,2,1},
{0,1,2},
{1,2,0},
{1,2,0},
};

//define face normals
const float face_normals[6][3] = {
{0.,0.,1.},
{0.,1.,0.},
{0.,-1.,0.},
{0.,0.,-1.},
{1.,0.,0.},
{-1.,0.,0.},
};

//define RGB for formal colors of face middles UBFDRL
const float cubergb[6][3] = {
{1.,1.,0.},
{0.,0.,1.},
{0.,1.,0.},
{1.,1.,1.},
{1.,.5,0.},
{1.,0.,0.},
};


const int facetext_rot[6][4]={
{0.,0.,0.,0.},
{-90.,1.,0.,0.},
{90.,1.,0.,0.},
{180.,1.,0.,0.},
{90.,0.,1.,0.},
{-90.,0.,1.,0.},
};

const int facetext_tr[6][3]={
{-10.,0.,0.},
{-10.,0.,0.},
{-10.,0.,0.},
{-10.,0.,0.},
{0.,10.,0.},
{0.,-10.,0.},
};




///////////////////////////////////////////////////////////////////////////////
//Rubik's group routines
///////////////////////////////////////////////////////////////////////////////


typedef int PERM[NP+1];

typedef struct 
{
PERM ep;
PERM vp;
} ELEMENT;


typedef int CYCLE[NP+1]; //element 0 is cycle length

typedef int cycle4[4];

typedef struct
{
cycle4 ec[2];
cycle4 vc[3];
} GENERATOR;


typedef struct
{
int generatornumber;
int generatorpower;
} GENERATORPOWER;


const GENERATOR Ug = {
{{1,2,3,4},{9,21,5,17}},{{1,2,3,4},{9,21,5,17},{12,24,8,20}}
};

const GENERATOR Fg = {
{{9,10,11,12},{13,22,3,20}},{{9,10,11,12},{13,22,4,20},{16,21,3,19}}
};

const GENERATOR Dg = {
{{13,14,15,16},{7,23,11,19}},{{13,14,15,16},{7,23,11,19},{6,22,10,18}}
};

const GENERATOR Bg = {
{{5,6,7,8},{1,24,15,18}},{{5,6,7,8},{2,24,15,18},{1,23,14,17}}
};

const GENERATOR Rg = {
{{17,18,19,20},{10,2,8,14}},{{17,18,19,20},{10,3,8,14},{9,2,7,13}}
};

const GENERATOR Lg = {
{{21,22,23,24},{4,12,16,6}},{{21,22,23,24},{12,16,6,1},{11,15,5,4}}
};


typedef int CYCLEPERM[NP][NP+1]; //number of cycles determined by zero-length cycles sentinel, cycle length is [][0] element

typedef struct
{
CYCLEPERM ec;
CYCLEPERM vc;
} CYCLEELEMENT;


void fprint_generator(FILE *f, const GENERATOR *g)
{
int ic,i;
for(ic=0; ic<2; ++ic) //perform edge cycles
        {
        for(i=0; i<4; ++i) fprintf(f,"%c%d%c",i==0?'(':' ',g->ec[ic][i],i==3?')':',');
        }
for(ic=0; ic<3; ++ic) //perform vertex cycles
        {
	for(i=0; i<4; ++i) fprintf(f,"%c%d%c",i==0?'(':' ',NP+g->vc[ic][i],i==3?')':',');
        }
}



void initperm(PERM p)
{
int i;
for(i=0; i<=NP; ++i) p[i]=i;
}

void initelement(ELEMENT *e)
{
initperm(e->ep);
initperm(e->vp);
}

void generator2element(const GENERATOR *g, ELEMENT *e)
{
initelement(e);
int ic,i;
for(ic=0; ic<2; ++ic) //perform edge cycles
	{
	for(i=0; i<4; ++i) e->ep[g->ec[ic][i]] = g->ec[ic][(i+1)%4];
	}
for(ic=0; ic<3; ++ic) //perform vertex cycles
        {
	for(i=0; i<4; ++i) e->vp[g->vc[ic][i]] = g->vc[ic][(i+1)%4];
        }
}




//returns pointer after closing bracket or NULL if no cycle found
//or input error
char *read1cycle(CYCLE c, const char *p)
{
c[0]=0;
if(*p==0) return NULL;
char *openbracket = strchr(p,'(');
if(!openbracket) return NULL;
char *closebracket = strchr(openbracket+1,')');
if(!closebracket) return NULL;
char *s = openbracket+1;
int r;
int nread=0;
int isvertex;
do	{
	int nchar;
	if(*s==',') ++s;
	r = sscanf(s,"%d%n",c+nread+1,&nchar);
	if(r==1)
		{
		++nread;
		s += nchar;
		if(c[nread]<1 || c[nread]> 2*NP) return NULL;
		if(nread==1) isvertex= c[1]>NP;
		else
			{
			if(c[nread]>NP && !isvertex || c[nread]<=NP && isvertex) return NULL;
			}
		}
	}
	while(r==1 && nread<NP && s<closebracket);
c[0]=nread;

/*
int i;
printf("read (");
for(i=1; i<=c[0]; ++i) printf("%d ",c[i]);
printf(")\n");
*/

return closebracket+1;
}

int ispermutation(PERM p)
{
PERM occured;
int i;
for(i=1; i<=NP; ++i) occured[i]=0;
for(i=1; i<=NP; ++i) occured[p[i]]++;
for(i=1; i<=NP; ++i) if(occured[i]!=1) return 0;
return 1;
}





void printcycleperm(const CYCLEPERM p, const int shift)
{
int i,j;
for(i=0; i<NP; ++i)
	{
	if(p[i][0])
		{
		printf("(");
		for(j=1;j<=p[i][0]; ++j) printf("%d%c",p[i][j]+shift,(j==p[i][0]?')':','));
		}
	else //terminated
		{
		if(i==0) printf("()");
		return;
		}
	}
}

void perm2cycles(CYCLEPERM c, const PERM p)
{
PERM used;
memset(used,0,sizeof(PERM));
int firstunused=1;
int currentcycle=0;
do
	{
	//find a cycle starting with first unused element
	int cyclelength=0;
	int next = firstunused;
	do
		{
		++cyclelength;
		used[next]=1;
		c[currentcycle][cyclelength] = next;
		next = p[next];
		}
	while(used[next]==0);
	if(cyclelength>1) //nontrivial cycle
		{
		c[currentcycle][0]=cyclelength;
		++currentcycle;
		}
	while(used[firstunused]) {++firstunused; if(firstunused>NP) break;} //find next unused element
	}
while(firstunused<=NP);
if(currentcycle<NP) c[currentcycle][0]=0;
}

void element2cycles(CYCLEELEMENT *c, const ELEMENT *e)
{
perm2cycles(c->ec,e->ep);
perm2cycles(c->vc,e->vp);
}


void printcycles(const CYCLEELEMENT *c)
{
printf(" edge  cycles: "); printcycleperm(c->ec,0);
printf("\nvertex cycles: "); printcycleperm(c->vc,NP);
printf("\n");
}

void printelement(const ELEMENT *e)
{
int j;
printf(" edge  permutation = ");
for(j=1; j<=NP; ++j) printf("%2d ",e->ep[j]);
printf("\n");
printf("vertex permutation = ");
for(j=1; j<=NP; ++j) printf("%2d ",e->vp[j]+NP);
printf("\n");

//print also in cycle form
CYCLEELEMENT c;
element2cycles(&c,e);
printcycles(&c);
}


void cycles2perm(PERM p, const CYCLEPERM c)
{
initperm(p);
int i,j; 
for(j=0; j<NP; ++j)
        {
        if(c[j][0] == 0) break;
        for(i=1; i<=c[j][0]; ++i) p[c[j][i]] = c[j][i%c[j][0]+1];
        }

if(!ispermutation(p)) printf("non-disjunct cycle input or internal error\n");
}

void cycles2element(ELEMENT *e, const CYCLEELEMENT *ce)
{
cycles2perm(e->vp,ce->vc);
cycles2perm(e->ep,ce->ec);
}


//do not check disjunct cycles here, but later in the conversion to explicit permutation
void scancycleelement(CYCLEELEMENT *ce, const char *p)
{
CYCLE c;
int i;
for(i=0; i<NP; ++i) ce->vc[i][0]=ce->ec[i][0]=0;

int ecount=0;
int vcount=0;
while(p=read1cycle(c,p)) 
	{
	//printf("cycle %d of length %d read\n",count,c[0]);
	if(c[0]==0) 
		{
		if(ecount==0 && vcount==0) return; 
		else  {printf("extra empty cycle\n"); exit(1);}
		} 
	//store a nontrivial cycle
	if(c[1]<=NP)
		{
		if(ecount>=NP) {printf("too many cycles\n"); return;}
		ce->ec[ecount][0] = c[0];
		for(i=1; i<=c[0]; ++i) ce->ec[ecount][i] = c[i];
		++ecount;
		}
	else
		{
		if(vcount>=NP) {printf("too many cycles\n"); return;}
		ce->vc[vcount][0] = c[0];
		for(i=1; i<=c[0]; ++i) ce->vc[vcount][i] = c[i]-NP;
		++vcount;
		}
	}
if(ecount+vcount==0) printf("empty or errorneous cycle input\n");
}


void scanpermelement(ELEMENT *e, const char *t)
{
const char *p = t;
int i,n;
for(i=1; i<=NP; ++i) 
	{
	if(1!=sscanf(p,"%d%n",&e->ep[i],&n))  {printf("input error\n"); return;}
	p+=n;
	}
for(i=1; i<=NP; ++i) 
        {
        if(1!=sscanf(p,"%d%n",&e->vp[i],&n))  {printf("input error\n"); return;}
	if(e->vp[i]<=NP) {printf("input error 2\n"); return;}
	e->vp[i] -= NP;
        p+=n;
        }
	

if(!ispermutation(e->ep) || !ispermutation(e->vp)) printf("input error - is not a permutation\n"); 
}


void scanelement(ELEMENT *e, const char *t)
{
char tmp[1024];
char *p = tmp;
while(*t)
	{
	if(*t == ','||*t == '\n') {*p++ = ' '; ++t;}
	else *p++ = *t++;
	}

if(strchr(tmp,'('))
	{
	CYCLEELEMENT c;
	scancycleelement(&c,tmp);
	//printf("test: ");printcycles(&c);
	cycles2element(e,&c);
	return;
	}
//scan plain permutation
scanpermelement(e,tmp);
}


//group operations

void perminv(PERM pinv, const PERM p)
{
int i;
for(i=1; i<=NP; ++i) pinv[p[i]]=i;
}

void permmult(PERM res, const PERM l, const PERM r)
{
int i;
for(i=1; i<=NP; ++i) res[i] = l[r[i]];
}

void groupinv(ELEMENT *einv, const ELEMENT *e)
{
perminv(einv->ep,e->ep);
perminv(einv->vp,e->vp);
}

void groupmult(ELEMENT *res, const ELEMENT *l,  const ELEMENT *r)
{
permmult(res->ep, l->ep, r->ep);
permmult(res->vp, l->vp, r->vp);
}

void groupconjugate(ELEMENT *res, const ELEMENT *l,  const ELEMENT *r) //lrl^-1
{
ELEMENT tmp,tmp2;
groupmult(&tmp,l,r);
groupinv(&tmp2,l);
groupmult(res,&tmp,&tmp2);
}

void groupcommutator(ELEMENT *res, const ELEMENT *l,  const ELEMENT *r) //l^-1r^-1lr
{
ELEMENT tmp,tmp2;
groupmult(&tmp,r,l);
groupinv(&tmp2,&tmp); // l^-1r^-1
groupmult(&tmp,l,r);
groupmult(res,&tmp2,&tmp);
}


int ispermidentity(const PERM p)
{
int i;
for(i=1; i<=NP; ++i) if(p[i]!=i) return 0;
return 1;
}

int ispermequal(const PERM p1,const PERM p2)
{
int i;
for(i=1; i<=NP; ++i) if(p1[i]!=p2[i]) return 0;
return 1;
}


int isidentity(const ELEMENT *e)
{
return ispermidentity(e->ep) && ispermidentity(e->vp);
}

int isequal(const ELEMENT *e1,const ELEMENT *e2)
{
return ispermequal(e1->ep,e2->ep) && ispermequal(e1->vp,e2->vp);
}



//compute formal colors[6][3][3] = 0(up) 1(back) 2(front) 3(down) 4(right) 5(left)
//by applying the edge and vertex permutation to initial colors color[i][x][y] = i
//also store which number is now at given position
//
void element2colors(const ELEMENT *state, int colors[6][3][3], int occupation[6][3][3])
{
int i;
ELEMENT invstate;
groupinv(&invstate,state);
for(i=0; i<6; ++i) //faces
	{
	int j,k;
	for(j=0; j<3; ++j) for(k=0; k<3; ++k) 
		{
		int n= cube_numbering[i][j][k];
		//find what is on that position by inverse permutation
		int c = -2;
		int p = -1;
		if(is_edge[j][k])
			{
			p = invstate.ep[n];
                        c = edgecolors[p];
                        }
		else if(is_vertex[j][k])
			{
			p = invstate.vp[n];
			c = vertexcolors[p];
			}
		else 
			{
			c = i; //face middles by definition
			p = 0;
			}
		if(c<0||p<0) {printf("internal error in color\n"); exit(1);}
		colors[i][j][k]=c;
		occupation[i][j][k]=p;
		}
	/*debug
	printf("colors of face %d (in internal, not geometrical order)\n",i);
	for(j=0; j<3; ++j) 
		{
		for(k=0; k<3; ++k) printf("%d ",colors[i][j][k]);
		printf("\n");
		}
	*/

	}
}


//compute order of an element
int orderof(const ELEMENT *e)
{
int i=1;
ELEMENT tmp= *e;
while(!isidentity(&tmp))
	{
	ELEMENT tmp2 = tmp;
	groupmult(&tmp,e,&tmp2);
	++i;
	}
return i;
}

static ELEMENT generators[6],invgenerators[6],sqgenerators[6];
static ELEMENT *generatoraddress[4]={NULL,generators,sqgenerators,invgenerators};


//apply sequence of generators (given by a string read from left)
void applygenerators(int inversed, ELEMENT *final, const int n, const GENERATORPOWER *s, const ELEMENT *initial)
{
int i;
ELEMENT tmp,tmp2;
tmp= *initial;
for(i=0; i<n; ++i)
	{
	int ii= inversed?n-i-1:i;
	tmp2=tmp;
	if(s[ii].generatorpower==0) tmp=tmp2;
	else
		{
		int p = s[ii].generatorpower;
		groupmult(&tmp,&generatoraddress[inversed?4-p:p][s[ii].generatornumber],&tmp2);
		}
	}
*final = tmp;
}


//string of name+power pairs like "RU2R-" applied left to right
int string2generatorsequence(GENERATORPOWER *g, int nmax, const char *s)
{
int n=0;
const char *p=s;
while(*p)
	{
	//fprintf(stderr,"parsing |%s|\n",p);
	switch(*p) {
	case 'u':
	case 'U':
		g[n].generatornumber=0;
		break;
        case 'b':
        case 'B':
                g[n].generatornumber=1;
                break;
        case 'f':
        case 'F':
                g[n].generatornumber=2;
                break;
        case 'd':
        case 'D':
                g[n].generatornumber=3;
                break;
        case 'r':
        case 'R':
                g[n].generatornumber=4;
                break;
        case 'l':
        case 'L':
                g[n].generatornumber=5;
                break;
	case ' ':
	case '\n':
	case 0:
		++p;
		continue;
	default: return -2;
	}
	if(p[1]==0||isalpha(p[1])) //just letter
		{
		g[n].generatorpower=1;
		++n;
        	if(n>=nmax) return -4;
        	++p;
		continue;
		}

	switch(p[1]){
	case '0':
		g[n].generatorpower=0;
		break;
	case '\n':
	case '1':
	case ' ':
	case '+':
		g[n].generatorpower=1;
                break;
	case '2':
		g[n].generatorpower=2;
                break;
	case '-':
	case '\'':
	case '3':
		g[n].generatorpower=3;
                break;
	default: return -3;
	}
	++n;
	if(n>=nmax) return -4;
	p+=2;
	}

return n;
}



void group_init(void)
{
char generatornames[6]="UFDBRL";
generator2element(&Ug,&generators[0]);
generator2element(&Bg,&generators[1]);
generator2element(&Fg,&generators[2]);
generator2element(&Dg,&generators[3]);
generator2element(&Rg,&generators[4]);
generator2element(&Lg,&generators[5]);

//check and print generators
#if 0
{
int i;
for(i=0; i<6; ++i)
	{
	if(!ispermutation(generators[i].ep) || !ispermutation(generators[i].vp)) {printf("internal error\n"); exit(1);}
	printf("Generator %c  = \n",generatornames[i]);
	printelement(generators+i);
	}
}
#endif

//prepare squares and inverses of the generators and perform consistency check

{
int i;
ELEMENT test,test2;
for(i=0; i<6; ++i)
        {
	groupinv(invgenerators+i,generators+i);
	groupmult(sqgenerators+i,generators+i,generators+i);
	groupmult(&test,invgenerators+i,generators+i);
	if(!isidentity(&test)) {printf("internal error\n"); exit(1);}
	groupmult(&test2,invgenerators+i,invgenerators+i);
	if(!isequal(&test2,sqgenerators+i)) {printf("internal error\n"); exit(1);}
	}
}

//check and print invgenerators
#if 0
{
int i;
for(i=0; i<6; ++i)
        {
        if(!ispermutation(invgenerators[i].ep) || !ispermutation(invgenerators[i].vp)) {printf("internal error\n"); exit(1);}
        printf("Inverse of Generator %c  = \n",generatornames[i]);
        printelement(invgenerators+i);
        }
for(i=0; i<6; ++i)
        {
        if(!ispermutation(sqgenerators[i].ep) || !ispermutation(sqgenerators[i].vp)) {printf("internal error\n"); exit(1);}
        printf("Square of Generator %c  = \n",generatornames[i]);
        printelement(sqgenerators+i);
        }

}
#endif

}

////////////////////////////////////////////////////////////////////////////////////////
//visualization routines
////////////////////////////////////////////////////////////////////////////////////////

const GLfloat cube_size=30., gap=5.;

GLfloat distance, aspect,window_w=1000,window_h=800;
GLint rot_x, rot_y;
const GLint rot_increment=3;


void init_visual(void)
{

//viewpoint
  rot_x = 0.0; 
  rot_y = 0.0; 
  distance = 50.;


  // init lighting
  GLfloat ambient_lighte[4]={0.2,0.2,0.2,1.0};
  GLfloat diffuse_light[4]={0.7,0.7,0.7,1.0};          
  GLfloat specular_light[4]={1.0, 1.0, 1.0, 1.0};   
  GLfloat light_position[4]={0.0, 50.0, 50.0, 1.0};

  // material brightness capacity
  GLfloat specularity[4]={1.0,1.0,1.0,1.0};
  GLint material_specularity = 60;

  // black background
  glClearColor(0.0f, 0.0f, 0.0f, 1.0f);

  // Gouraud colorization model
  glShadeModel(GL_SMOOTH);

  // material reflectability
  glMaterialfv(GL_FRONT,GL_SPECULAR, specularity);
  // brightness concentration
  glMateriali(GL_FRONT,GL_SHININESS,material_specularity);

  // activate ambient light
  glLightModelfv(GL_LIGHT_MODEL_AMBIENT, ambient_lighte);

  // define light parameters
  glLightfv(GL_LIGHT0, GL_AMBIENT, ambient_lighte);
  glLightfv(GL_LIGHT0, GL_DIFFUSE, diffuse_light );
  glLightfv(GL_LIGHT0, GL_SPECULAR, specular_light );
  glLightfv(GL_LIGHT0, GL_POSITION, light_position );

  // enable changing material color
  glEnable(GL_COLOR_MATERIAL);
  // enable lighting
  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);
  // enable depth buffering
  glEnable(GL_DEPTH_TEST);
}


int cube_colors[6][3][3];
int cube_occupations[6][3][3];
int plottext;


void draw_callback(void)
{
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glLoadIdentity();
gluLookAt(0,80,200, 0,0,0, 1,1,1);

//rotate the cube as whole
glRotatef(rot_x, 1.0, 0.0, 0.0);
glRotatef(rot_y, 0.0, 1.0, 0.0);

//draw the cube
int i;
for(i=0; i<6; ++i) //faces
	{
	//prepare coordinates	
	int coord1=face_planes_and_normal[i][0];
	int coord2=face_planes_and_normal[i][1];
	int normal=face_planes_and_normal[i][2];
	GLfloat fnormal=face_normals[i][normal];
	int j,k;
	GLfloat xyz[2][2][3];
	for(j=0; j<2; ++j) for(k=0; k<2; ++k) //coordinates of one cubicle
		{
		xyz[j][k][coord1] = (1.*j-.5)*cube_size;
		xyz[j][k][coord2] = (1.*k-.5)*cube_size;
		xyz[j][k][normal] =  (1.5*cube_size+gap)*fnormal;
		}

	//draw all cubicles
	int ic,jc;
	for(ic=0; ic<3; ++ic) for(jc=0; jc<3; ++jc) //loop over cubicles
		{
		int c; //formal color
		c= cube_colors[i][ic][jc];
		if(c<0||c>5) {printf("wrong color\n"); exit(1);}
		glPushMatrix(); 
		glColor3f(cubergb[c][0],cubergb[c][1],cubergb[c][2]);
		
		//translate the little cubicle
		GLfloat translation[3];
		translation[coord1] = (ic-1)*(cube_size+gap);
		translation[coord2] = (jc-1)*(cube_size+gap);
		translation[normal] = 0.;
		glTranslatef(translation[0],translation[1],translation[2]); 

		//draw square
		glBegin(GL_QUADS);
		glNormal3f(face_normals[i][0],face_normals[i][1],face_normals[i][2]);
		glVertex3f(xyz[0][0][0],xyz[0][0][1],xyz[0][0][2]);
		glVertex3f(xyz[0][1][0],xyz[0][1][1],xyz[0][1][2]);
		glVertex3f(xyz[1][1][0],xyz[1][1][1],xyz[1][1][2]);
		glVertex3f(xyz[1][0][0],xyz[1][0][1],xyz[1][0][2]);
		glEnd();
		glPopMatrix();
		
		//text labels
		if(plottext)
		{
		glPushMatrix();
		translation[coord1] = (ic-1)*(1.02*cube_size+gap)+facetext_tr[i][0];
                translation[coord2] = (jc-1)*(1.02*cube_size+gap)+facetext_tr[i][1];
                translation[normal] =  (1.52*cube_size+gap)*fnormal+facetext_tr[i][2];
                glTranslatef(translation[0],translation[1],translation[2]);
		char label[16];
		int shift = is_vertex[ic][jc]?NP:0;
		if(ic==1&&jc==1) sprintf(label,"  %c  ",faces[i]);
		else sprintf(label,"%2d@%2d",cube_occupations[i][ic][jc]+shift,cube_numbering[i][ic][jc]+shift);
		glScalef(0.04,0.04,0.04);
		glLineWidth(4);
		glColor3f(1.-cubergb[c][0],1.-cubergb[c][1],1.-cubergb[c][2]);
		glRotatef(facetext_rot[i][0],facetext_rot[i][1],facetext_rot[i][2],facetext_rot[i][3]);
		glutStrokeString(GLUT_STROKE_MONO_ROMAN,label);
		glPopMatrix();
		}//plottext
		}//for ic,jc

	//draw four black gaps
	glColor3f(0.,0.,0.);
	GLfloat xyz2[3];
	xyz2[normal] =  (1.5*cube_size+gap)*fnormal;

	glBegin(GL_QUADS);
        glNormal3f(face_normals[i][0],face_normals[i][1],face_normals[i][2]);
	xyz2[coord1] = -1.5*cube_size-gap;
	xyz2[coord2] = 0.5*cube_size+gap;
	glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
	xyz2[coord1] = -xyz2[coord1];
	glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
	xyz2[coord2] = 0.5*cube_size;
	glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
	xyz2[coord1] = -xyz2[coord1];
	glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
	glEnd();

	glBegin(GL_QUADS);
        glNormal3f(face_normals[i][0],face_normals[i][1],face_normals[i][2]);
        xyz2[coord1] = -1.5*cube_size-gap;
        xyz2[coord2] = -0.5*cube_size-gap;
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord1] = -xyz2[coord1];
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord2] = -0.5*cube_size;
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord1] = -xyz2[coord1];
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        glEnd();

	glBegin(GL_QUADS);
        glNormal3f(face_normals[i][0],face_normals[i][1],face_normals[i][2]);
        xyz2[coord2] = -1.5*cube_size-gap;
        xyz2[coord1] = 0.5*cube_size+gap;
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord2] = -xyz2[coord2];
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord1] = 0.5*cube_size;
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord2] = -xyz2[coord2];
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        glEnd();

        glBegin(GL_QUADS);
        glNormal3f(face_normals[i][0],face_normals[i][1],face_normals[i][2]);
        xyz2[coord2] = -1.5*cube_size-gap;
        xyz2[coord1] = -0.5*cube_size-gap;
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord2] = -xyz2[coord2];
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord1] = -0.5*cube_size;
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        xyz2[coord2] = -xyz2[coord2];
        glVertex3f(xyz2[0],xyz2[1],xyz2[2]);
        glEnd();



	}


glutSwapBuffers();
}

int mousex0,mousey0;
int mouse_left_pressed=0;
void mouse_callback(int button, int state, int x, int y)
{
if (button == GLUT_LEFT_BUTTON)
	{
	if (state == GLUT_DOWN) mouse_left_pressed=1;
	else mouse_left_pressed=0;
	mousex0=x;
	mousey0=y;
	}

glutPostRedisplay();
}



void motion_callback(int x, int y)
{
if(mouse_left_pressed)
	{
	//relative to x0, set rotation
	rot_x = (x-mousex0) % 360;
	rot_y = -(y-mousey0) % 360;
	}


glutPostRedisplay();
}

void keyboard_callback(unsigned char key, int x, int y)
{
switch(key) {
	case 'q':
		exit(0);
		break;
	default:
		break;
}
glutPostRedisplay();
}


void reshape_callback(GLsizei w, GLsizei h)
{
if ( h == 0 ) h = 1;
if ( w == 0 ) w = 1;
window_w=w; window_h=h;
glViewport(0, 0, w, h);
aspect = (GLfloat)w/(GLfloat)h;
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluPerspective(distance,aspect,0.4,500);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluLookAt(0,80,200, 0,0,0, 1,1,1);
}



void visualize(int argc, char **argv)
{
int pid=fork();
if(pid) return; //parent

  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutInitWindowSize(window_w,window_h);
  aspect=window_w/window_h;
  gluPerspective(distance,aspect,0.4,500);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  glutCreateWindow("Rubik 3D");
  glutDisplayFunc(draw_callback);
  glutReshapeFunc(reshape_callback);
  glutMouseFunc(mouse_callback);
  glutMotionFunc(motion_callback);
  glutKeyboardFunc(keyboard_callback);
  init_visual();
  glutMainLoop();
}


////////////////////////////////////////////////////////////////////////////////
//GAP interface routines
///////////////////////////////////////////////////////////////////////////////
FILE *from_gap, *to_gap;
int GAP_available=0;

void pipe_handler(int i)
{
GAP_available=0;
}

//try to open the GAP program for advanced operations
void GAP_init(void)
{
int pid;
int pipes[2][2];

if(pipe(pipes[0])) perror("pipe error");
if(pipe(pipes[1])) perror("pipe error");

GAP_available=1;
signal(SIGPIPE,pipe_handler);

pid=fork();

if(pid==0) //child
  {
  close(pipes[0][0]);
  close(pipes[1][1]);

  int c1=dup2(pipes[0][1],STDOUT_FILENO);
  if(c1<0) perror("dup2 error");
  close(pipes[0][1]);
  int c2=dup2(pipes[1][0],STDIN_FILENO);
  if(c2<0) perror("dup2 error");
  close(pipes[1][0]);

   setlinebuf(stdin);
   setlinebuf(stdout);


  char *(gap_argv[64]) = {"-b","-q","-E","-T",NULL};
  if(execvp("gap",gap_argv)) {printf("GAP NOT FOUND\n"); exit(1);};
  exit(0);
  }

//parent
close(pipes[0][1]);
close(pipes[1][0]);


from_gap = fdopen(pipes[0][0],"r");
to_gap = fdopen(pipes[1][1],"w");
if(from_gap==NULL || to_gap==NULL) {GAP_available=0; return;}

//line buffering to avoid need to fflush()
setlinebuf(from_gap);
setlinebuf(to_gap);

return;
}


void GAP_setup(void)
{
if(!GAP_available) return;
char *gapout; char buf[1024];
int i;
fprintf(to_gap,"cube := Group( ");
fprint_generator(to_gap,&Ug);fprintf(to_gap,",");
fprint_generator(to_gap,&Bg);fprintf(to_gap,",");
fprint_generator(to_gap,&Fg);fprintf(to_gap,",");
fprint_generator(to_gap,&Dg);fprintf(to_gap,",");
fprint_generator(to_gap,&Rg);fprintf(to_gap,",");
fprint_generator(to_gap,&Lg);
fprintf(to_gap,");\n");

//pipe signal might have happened if GAP not available
if(!GAP_available) return;

gapout = fgets(buf,1024,from_gap); //if(gapout) fputs(buf,stdout);
fprintf(to_gap,"f := FreeGroup(\"U\",\"B\",\"F\",\"D\",\"R\",\"L\");\n");
gapout = fgets(buf,1024,from_gap); //if(gapout) fputs(buf,stdout);
fprintf(to_gap,"hom := GroupHomomorphismByImages( f, cube, GeneratorsOfGroup(f),GeneratorsOfGroup(cube) );\n");
gapout = fgets(buf,1024,from_gap); //if(gapout) fputs(buf,stdout);
for(i=1; i<=6; ++i) {gapout = fgets(buf,1024,from_gap); /*if(gapout) fputs(buf,stdout);*/}

}


void GAP_shell(char *l)
{
if(!GAP_available) return;
fprintf(to_gap,"%s\n",l);
struct timespec delay = {0,600000000};
nanosleep(&delay,&delay);

//set non-blocking I/O
int flags;
flags = fcntl(fileno(from_gap), F_GETFL, 0);
flags |= O_NONBLOCK;
fcntl(fileno(from_gap), F_SETFL, flags);

//read all stuff from GAP
char *gapout;
do
  {
  char buf[1024];
  gapout = fgets(buf,1024,from_gap);
  if(gapout)
        {
        fputs(buf,stdout);
        }
  }
while(gapout);

//set blocking I/O again
flags = fcntl(fileno(from_gap), F_GETFL, 0);
flags &= ~O_NONBLOCK;
fcntl(fileno(from_gap), F_SETFL, flags);
}


int GAP_ingroup(ELEMENT *e)
{
if(!GAP_available) return 0;
int j;
fprintf(to_gap,"AsPermutation(Transformation([");
for(j=1; j<=NP; ++j) fprintf(to_gap,"%d,",e->ep[j]);
for(j=1; j<=NP; ++j) fprintf(to_gap,"%d%c",e->vp[j]+NP,j<NP?',':' ');
fprintf(to_gap,"]));\n last in cube;\n");

struct timespec delay = {0,500000000};
nanosleep(&delay,&delay);

//read from GAP
char *gapout, buf[1024];
int found= -1;
while(found== -1)
        {
	gapout = fgets(buf,1024,from_gap);
	if(gapout && strstr(buf,"true")) found=1;
	if(gapout && strstr(buf,"false")) found=0;
        }
return found;
}


void GAP_random(ELEMENT *e)
{
char *gapout;
char buf[1024],tmp[2048];
if(!GAP_available) {printf("GAP not available\n"); return;}

int flags;
flags = fcntl(fileno(from_gap), F_GETFL, 0);
flags &= ~O_NONBLOCK;
fcntl(fileno(from_gap), F_SETFL, flags);

fprintf(to_gap,"Random(cube);\n");
gapout = fgets(buf,1024,from_gap); 
if(gapout) strcpy(tmp,buf);

{struct timespec delay = {0,300000000}; nanosleep(&delay,&delay);}
//set non-blocking I/O
flags = fcntl(fileno(from_gap), F_GETFL, 0);
flags |= O_NONBLOCK;
fcntl(fileno(from_gap), F_SETFL, flags);
//read all remaining stuff from GAP
do
  {
  gapout = fgets(buf,1024,from_gap);
  if(gapout) strcat(tmp,buf);
  }
while(gapout);

//printf("random = %s\n",tmp);
scanelement(e,tmp);
}




char *GAP_decompose(char *string, ELEMENT *e)
{
string[0]=0;
char *gapout;
char buf[1024];
if(!GAP_available) return string;
int j;
fprintf(to_gap,"AsPermutation(Transformation([");
for(j=1; j<=NP; ++j) fprintf(to_gap,"%d,",e->ep[j]);
for(j=1; j<=NP; ++j) fprintf(to_gap,"%d%c",e->vp[j]+NP,j<NP?',':' ');
fprintf(to_gap,"]));\n");

gapout = fgets(buf,1024,from_gap); //if(gapout) fputs(buf,stdout); //read one line

{struct timespec delay = {0,300000000}; nanosleep(&delay,&delay);}
//set non-blocking I/O
int flags;
flags = fcntl(fileno(from_gap), F_GETFL, 0);
flags |= O_NONBLOCK;
fcntl(fileno(from_gap), F_SETFL, flags);
//read all remaining stuff from GAP
do
  {
  gapout = fgets(buf,1024,from_gap);
  //if(gapout) fputs(buf,stdout);
  }
while(gapout);


fprintf(to_gap,"PreImagesRepresentative(hom,last);\n");
{struct timespec delay = {0,900000000}; nanosleep(&delay,&delay);}

flags = fcntl(fileno(from_gap), F_GETFL, 0);
flags &= ~O_NONBLOCK;
fcntl(fileno(from_gap), F_SETFL, flags);
gapout = fgets(buf,1024,from_gap); //if(gapout) fputs(buf,stdout); //read one line
if(gapout && strstr(buf,"identity")) return string;
if(gapout) strcat(string,buf);

flags = fcntl(fileno(from_gap), F_GETFL, 0);
flags |= O_NONBLOCK;
fcntl(fileno(from_gap), F_SETFL, flags);
//read remainder
do
  {
  gapout = fgets(buf,1024,from_gap);
  //if(gapout) fputs(buf,stdout);
  if(gapout) strcat(string,buf);
  }
while(gapout);

return string;
}

//assumes no nested parenthesis and single digit exponents and syntactically coorect output from GAP for simplicity
char *reformat_GAP_expression(char *out, char *in)
{

char *p,*q,*r,*s,*t;
char tmp[2048];
p=tmp; q=in;

while(*q) //get rid of some redundant signs
        {
        if(*q=='\\' ||*q=='\n' || *q=='*' || *q=='^') ++q;
	else *p++ = *q++;
	}
*p=0;


//perform simple parentheses processing
char tmp2[2048];
p=tmp2; q=tmp; r=NULL;
while(*q)
	{
	if(*q=='(')
		{
		++q;
		r=p;
		}
	else if(*q==')')
		{
		++q;
		int nrep= *q-'0' - 1;
		if(nrep<0 || r==NULL||*q == '-') {printf("internal error in reformat_GAP_expression\n"); exit(1);}
		++q;
		int i;
		char *s= p;
		for(i=0; i<nrep; ++i) //copy the part nrep times to output
			{
			char *t;
			for(t=r; t<s; ++t) *p++ = *t;
			}
		r=NULL;
		}
	else *p++ = *q++;
	}
*p=0;


//canonicalize exponents
p=out; q=tmp2;
while(*q)
	{
	if(isalpha(*q)) *p++ = *q++;
	else if(isdigit(*q)) {if(*q=='3') *p++ = '-'; else if(*q=='2') *p++ = '2'; else {printf("internal error 3 in reformat_GAP_expression\n"); exit(1);} ++q;}
	else if(*q == '-') {++q; if(*q=='3') {;} else if(*q=='2') *p++ = '2'; else if(*q=='1') *p++ = '-'; else {printf("internal error 4 in reformat_GAP_expression\n"); exit(1);} ++q;} 
	else {printf("internal error 2 in reformat_GAP_expression\n"); exit(1);}
	}
*p=0;

return out;
}

#define MAXSEQUENCE 256

ELEMENT storage[10];
int storageOK(int n)
{
if(storage[n].ep[1]==0) {printf("Memory %d undefined\n",n); return 0;}
return 1;
}

int main(int argc, char **argv)
{
//for readline
rl_bind_key('\t',rl_insert);
rl_variable_bind("meta-flag","On");
rl_variable_bind("convert-meta","Off");
rl_variable_bind("output-meta","On");
rl_instream=stdin;
rl_outstream=stderr;

group_init();

GAP_init();
GAP_setup();
if(!GAP_available) printf("GAP not found, advanced functionality not available\n");

ELEMENT currentstate;
memset(storage,0,sizeof(storage));
initelement(&currentstate);

while(1)
	{
	char *line;
	line = readline("|rubik> ");
	if(!line) goto end;
	if(*line == 0 || *line == '\n') continue;
	else add_history(line);
	char *l=line;
	while(isblank(*l)) ++l;
	if(*l) switch(*l) {
		case 'g':
			printf("%d\n",GAP_ingroup(&currentstate));
			break;
		case 's':
			{char buf[2048]; char buf2[1024];
			GAP_decompose(buf, &currentstate);
			printf("%s\n",reformat_GAP_expression(buf2,buf));
			}
			break;
		case '!':
			GAP_shell(l+1);
			break;
		case 'w':
			if(isdigit(l[1])) storage[l[1]-'0'] = currentstate;
			break;
		case 'r':
			if(isdigit(l[1])) 
				{
				if(storageOK(l[1]-'0')) currentstate = storage[l[1]-'0'];
				}
                        break;
		case 'I':
			printf("%d\n",isidentity(&currentstate));
			break;
		case 'i':
			if(isdigit(l[1])) groupinv(&currentstate, storage + l[1]-'0');
                        break;
		case 'e':
			if(isdigit(l[1]) && isdigit(l[2]) && storageOK(l[1]-'0') && storageOK(l[2]-'0')) printf("%d\n",isequal(storage + l[1]-'0', storage + l[2]-'0'));
			break;
		case 'm':
			if(isdigit(l[1]) && isdigit(l[2]) && storageOK(l[1]-'0') && storageOK(l[2]-'0')) groupmult(&currentstate, storage + l[1]-'0', storage + l[2]-'0');
			break;
		case 'c':
			if(isdigit(l[1]) && isdigit(l[2]) && storageOK(l[1]-'0') && storageOK(l[2]-'0'))  groupconjugate(&currentstate, storage + l[1]-'0', storage + l[2]-'0');
                        break;
		case 'C':
			if(isdigit(l[1]) && isdigit(l[2]) && storageOK(l[1]-'0') && storageOK(l[2]-'0'))  groupcommutator(&currentstate, storage + l[1]-'0', storage + l[2]-'0');
			break;
		case 'X':
		case 'x':
			{
			ELEMENT tmpstate;
			GENERATORPOWER seqcode[MAXSEQUENCE];
			int nseq=string2generatorsequence(seqcode,MAXSEQUENCE,l+1);
			if(nseq<=0) printf("Empty or illegal operation %d: %s\n",nseq,l+1);
			//int i; for(i=0; i<nseq; ++i) printf("Sequence %d %d\n",seqcode[i].generatornumber,seqcode[i].generatorpower);
			applygenerators(*l=='X',&tmpstate,nseq,seqcode,&currentstate);
			currentstate=tmpstate;
			//printelement(&currentstate);
			}
			break;
		case 'v':
			if(l[1]=='0') plottext=0; else plottext=1;
			element2colors(&currentstate,cube_colors,cube_occupations);
			visualize(argc,argv);
			break;
		case '1':
			printf("state reset\n");
			initelement(&currentstate);
			break;
		case 'R':
			GAP_random(&currentstate);
			printf("randomized\n");
			break;
		case 'P':
			scanelement(&currentstate,l+1);
			break;
		case 'p':
			printelement(&currentstate);
			break;
		case 'o':
			{
			printf("%d\n",orderof(&currentstate));
			}
			break;
		case 'q':
			goto end;
			break;
		case '\n':
			break;
		case 'h':
		case 'H':
		case '?':
			printf("Commands:\n");
			printf("1 ... initialize state to solved cube\n");
			printf("P ... input permutation (1-48) straight or in parenthesized cycle form\n");
			printf("R ... generate random permutation from the cube group\n");
			printf("p ... print current permutation\n");
			printf("v ... visualize current state (v0 without permutation labels)\n");
			printf("x sequence ... perform sequence of operations, e.g. x r2l-ubd2f\n");
			printf("X sequence ... perform inverse of the sequence\n");
			printf("o ... find the order (period) of the current element\n");
			printf("w1 ... store state to memory 1\n");
			printf("r9 ... read state from memory 9\n");
			printf("i3 ... compute inverse of memory 3\n");
			printf("I ... is identity?\n");
			printf("c01 ... conjugate memory0 memory1 memory0^-1\n");
			printf("C01 ... commutation of memory0 memory1\n");
			printf("m01 ... multiply memory0 * memory1\n");
			printf("! ... escape to GAP shell\n");
			printf("g ... check if the current permutation is in the cube group\n");
			printf("s ... solve the current position - decompose to generator sequence\n");
			printf("h ... help\n");
			printf("q ... quit\n");
			break;
		default:
			printf("unknown command\n");
			break;
	}//switch
	if(line) free(line);
	}//while

end:;
killpg(0,SIGTERM);
int status;
wait(&status);
return 0;
}