#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <complex.h>
#include <string.h>

#include "../numeric/SU21-numeric.h"
#include "../numeric/SU21-numeric-base.c"

#define GOODPREFIX "C1/C1-5"

// \sqrt{5}
#define sq5 (2.236067977499789696409173668731)
// The golden ration
#define GR (1.618033988749894848204586834366)
// After one round of conjugation, the generators act unitarily with
// respect to the form $\diag(1,1,-\GR)=\diag(1,1,-R_0^2)$.

// R0 =$\sqrt(\GR))$
#define R0 (1.272019649514068964252422461737)

// Powers of $Z=\zeta_5=e^{2\pi i/5}$
#define Z  ( 0.309016994374947424102293417183+0.951056516295153572116439333379*I)
#define Z2 (-0.809016994374947424102293417183+0.587785252292473129168705954639*I)
#define Z3 (-0.809016994374947424102293417183-0.587785252292473129168705954639*I)
#define Z4 ( 0.309016994374947424102293417183-0.951056516295153572116439333379*I)

typedef struct {
  char Nam[MaxGenNameLength+1]; // Generator name
  char InvNam[MaxGenNameLength+1]; // Inverse generator name
  int Len; // ``Length'' of Generator
  Mtx_t Mtx; // matrix for generator
} Gen_t;

#define MaxNoGens 1
#define ActualNoGens 1
#define KSiz 600

void SetupGens(Elt_t GenElts[],int *NoGenElts_p,
	       Mtx_t FF0,Mtx_t ConjMtx,Mtx_t ConjMtxInv);
void SetupK(Elt_t KElts[],
	       Mtx_t FF0,Mtx_t ConjMtx,Mtx_t ConjMtxInv);

int main(int argc,char *(argv[])) {
  // The generators as group elements
  Elt_t GenElts[2*ActualNoGens],KElts[KSiz];

  char *Prefix_p,FileNam[MAX_FILENAME_LEN]="../";
  int Ix,GenIx,GenEltIx,NoGenElts;
  size_t NamLen;

  // The form with respect to which the original, given matrices are
  // unitary:
  Mtx_t FF0=
    {{ -sq5,  Z4-Z, 0},
     { Z-Z4,  -sq5, 0},
     {    0,     0, 1}};

  // First matrix by which to conjugate them; this makes them unitary with
  // respect to $\diag{1,1,-\GR}=\diag(1,1,-R_0^2)$
  Mtx_t ConjMtx1=
    {{ -Z+Z4, -Z2-Z3, 0 },
     { 0,         -1, 0 },
     { 0,          0, 1 }};

  // Second matrix by which to conjugate them to make them unitary with
  // respect to $\diag{1,1,-1}$
  Mtx_t ConjMtx2=
    {{ 1, 0,  0},
     { 0, 1,  0},
     { 0, 0, R0}};

  // Full matrix by which to conjugate them to make them unitary with
  // respect to $\diag{1,1,-1}$
  Mtx_t ConjMtx21,ConjMtx21Inv;

  FILE *GensFile;
  FILE *KFile;

  // Direction in which to shift the origin:
  //Pt_t Dir={0.5237+0.7208*I,0.1403+0.4318*I};
  Pt_t Dir={0.1855+0.1348*I,0.7862+0.5712*I};
  //Pt_t Dir={0.7862,0.4859};
  FLOAT_t sShift=0.800;   // How far to shift it (Euclidean distance)
  Pt_t W;

  // Check that the prefix is right
  Prefix_p=getenv("PREFIX");
  assert(Prefix_p != NULL);
  assert(strncmp(Prefix_p,GOODPREFIX,50) == 0);

  // Open the output files
  {
    size_t len=strlen(Prefix_p);
    assert(3+len+5<MAX_FILENAME_LEN);
    strcpy(&FileNam[3],Prefix_p);
    strcpy(&FileNam[3]+len,"-Gens");
    printf("%s\n",FileNam);
    GensFile=fopen(FileNam,"w");
    assert(GensFile != NULL);
    strcpy(&FileNam[3]+len,"-K");
    printf("%s\n",FileNam);
    KFile=fopen(FileNam,"w");
    assert(KFile != NULL);
  }

  // Calculate ConjMtx21 and its inverse
  MultMtx(ConjMtx21,ConjMtx2,ConjMtx1);
  InvMtx(ConjMtx21Inv,ConjMtx21);

  // Prepare the generators and write the generator file
  SetupGens(GenElts,&NoGenElts,FF0,ConjMtx21,ConjMtx21Inv);
  fprintf(GensFile,"NoGenElts=%i\n",NoGenElts);
  for(GenEltIx=0;GenEltIx<NoGenElts;GenEltIx++)
    OutputElt(GensFile,GenElts[GenEltIx]);

  // Prepare K (the finite group) and write the K file
  SetupK(KElts,FF0,ConjMtx21,ConjMtx21(Inv);
  fprintf(KFile,"KSiz=%i\n",KSiz);
  for(GenEltIx=0;GenEltIx<KSiz;GenEltIx++)
    OutputElt(KFile,KElts[GenEltIx]);
}

// Prepare the generating elements
//
void SetupGens(Elt_t GenElts[],int *NoGenElts_p,
	       Mtx_t FF0,Mtx_t ConjMtx,Mtx_t ConjMtxInv) {

  // The generators
  Gen_t Gens[MaxNoGens]={
    {"A","a",100,
     {
       {          -1,             0,       0},
       { -2*Z2-2*Z-1,      -Z3-Z2+1,  -(Z+1)},
       {   -Z2-2*Z-2, Z*(-Z2-3*Z-1), Z3+Z2-1}
     }
    }
  };

  int GenIx,GenEltIx;
  size_t NamLen;
  Mtx_t MtxInv,tmpMtx;

  GenEltIx=0;
  for(GenIx=0;GenIx<ActualNoGens;GenIx++) {
    // Check that the matrix is unitary with respect to FF0
    printf("Checking GenIx=%i\n",GenIx);
    // OutputMtx(stdout,Gens[GenIx].Mtx);
    assert(CheckPU(Gens[GenIx].Mtx,FF0,epsBig));
    // Calculate the conjugated matrix
    ConjugateMtx(GenElts[GenEltIx].Mtx,
		 Gens[GenIx].Mtx,ConjMtx,ConjMtxInv);
    // and check that that is unitary relative to FF1
    assert(CheckPU(GenElts[GenEltIx].Mtx,FF1,epsBigger));
    InvMtx(MtxInv,GenElts[GenEltIx].Mtx);
    // Set up the generator name
    NamLen=strlen(Gens[GenIx].Nam);
    assert(NamLen>0 && NamLen<MaxGenNameLength);
    strcpy(GenElts[GenEltIx].Word,Gens[GenIx].Nam);
    // and ``length''
    GenElts[GenEltIx].Len=Gens[GenIx].Len;
    // and KIx and KStart
    GenElts[GenEltIx].KIx = -1;
    GenElts[GenEltIx].KStart = 0;
    
    FixUpElt(&GenElts[GenEltIx]);
    GenEltIx++;

    // Now do the inverse generator.  If the inverse name is empty, we
    // skip this.
    NamLen=strlen(Gens[GenIx].InvNam);
    if(NamLen>0) {
      // The matrix
      memcpy(&GenElts[GenEltIx].Mtx,MtxInv,sizeof(Mtx_t));
      assert(CheckPU(GenElts[GenEltIx].Mtx,FF1,epsBigger));
      // the inverse generator name
      assert(NamLen<MaxGenNameLength);
      strcpy(GenElts[GenEltIx].Word,Gens[GenIx].InvNam);
      // the length
      GenElts[GenEltIx].Len=Gens[GenIx].Len;
      // and KIx and KStart
      GenElts[GenEltIx].KIx = -1;
      GenElts[GenEltIx].KStart = 0;

      FixUpElt(&GenElts[GenEltIx]);
      GenEltIx++;
    }
  }
  *NoGenElts_p=GenEltIx;
}

// Set up the 600-element finite group, K
//
void SetupK(Elt_t KElts[],
	    Mtx_t FF0,Mtx_t ConjMtx,Mtx_t ConjMtxInv) {

  int Ix1,Ix2,Ix3;
  Elt_t *KElt_p;
  Mtx_t tmpMtx,KMtx;

  // Matrix M --- not used directly
  // Mtx_t M={{ -1, Z3, 0},
  //          {-Z2,  0, 0},
  // 	      {  0,  0, 1}};
  // Matrix N and its powers
  Mtx_t N={{1, Z2, 0},
	   {0, Z4, 0},
	   {0,  0, Z}};
  Mtx_t NPow[5];
  char NPowNam[5][4]={"","N","N^2","N^3","N^4"};
  // Matrix J and its powers; J = (M*N^-1)^2
  Mtx_t J={{-Z,0,0},{0,-Z,0},{0,0,Z3}},JPow[10];
  char JPowNam[10][4]=
    {"","J","J^2","J^3","J^4","J^5","J^6","J^7","J^8","J^9"};
  // Twelve coset representatives for
  //    < N > \ K_{600} / < J >
  Mtx_t ListOf12[12]={
    {{1,0,0},{0,1,0},{0,0,1}},
    {{-1,Z3,0},{-Z2,0,0},{0,0,1}},
    {{-1,0,0},{-Z2,-Z4,0},{0,0,Z}},
    {{-1,-Z2,0},{-Z2,-Z3-Z4,0},{0,0,Z2}},
    {{-1,-Z-Z2,0},{-Z2,-Z2-Z3-Z4,0},{0,0,Z3}},
    {{-1,Z3+Z4,0},{-Z2,1,0},{0,0,Z4}},
    {{1,-Z3,0},{Z+Z2,-1,0},{0,0,Z}},
    {{1,0,0},{Z+Z2,Z3,0},{0,0,Z2}},
    {{1,Z2,0},{Z+Z2,Z2+Z3+Z4,0},{0,0,Z3}},
    {{1,Z+Z2,0},{Z+Z2,Z+Z2+2*Z3+Z4,0},{0,0,Z4}},
    {{1,-Z3-Z4,0},{Z+Z2,Z2+Z3,0},{0,0,1}},
    {{Z3+Z4,-Z-Z2-Z3,0},{-Z-2*Z2-Z3-Z4,-Z3-Z4,0},{0,0,Z4}}
  };

  char ListOf12Nam[12][10]={
    "","M","MN","MN^2","MN^3","MN^4",
    "MNM","MNMN","MNMN^2","MNMN^3","MNMN^4","MNMN^3M^2"
  };

  // Set up powers of J
  memcpy(JPow[0],Id3,sizeof(Mtx_t));
  for(Ix1=0;Ix1<9;Ix1++)
    MultMtx(JPow[Ix1+1],J,JPow[Ix1]);
  // Set up powers of N
  memcpy(NPow[0],Id3,sizeof(Mtx_t));
  for(Ix3=0;Ix3<4;Ix3++)
    MultMtx(NPow[Ix3+1],N,NPow[Ix3]);

  // Set up elements of $K_{600}$
  KElt_p=KElts;
  for(Ix1=0;Ix1<10;Ix1++)
    for(Ix2=0;Ix2<12;Ix2++) {
      MultMtx(tmpMtx,ListOf12[Ix2],JPow[Ix1]);
      for(Ix3=0;Ix3<5;Ix3++) {
	// Set up the matrix
	MultMtx(KMtx,NPow[Ix3],tmpMtx);
	assert(CheckUnitary(KMtx,FF0));
	ConjugateMtx(KElt_p->Mtx,KMtx,ConjMtx,ConjMtxInv);
	assert(CheckPU(KElt_p->Mtx,FF1,epsSmall));
	// the generator name
	sprintf(KElt_p->Word,"%s%s%s",
		NPowNam[Ix3],ListOf12Nam[Ix2],JPowNam[Ix1]);
	// the ``length''
	if(Ix1 == 0 && Ix2 ==0 && Ix3 == 0)
	  KElt_p->Len = 0;
	else
	  KElt_p->Len = 2;
	// and KIx and KStart
	KElt_p->KIx = KElt_p->KStart = (KElt_p-KElts);

	FixUpElt(KElt_p);
	KElt_p++;
      }
    }
}