#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 "C18p3/C18p3-2"

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;

// $\sqrt{6}$ 
#define sq6 (2.449489742783178098197284074706)
// $\sqrt{-3}$ 
#define sqm3 (1.732050807568877293527446341506*I)
// $e^{2\pi i}/9$, $e^{4\pi i}/9$, $e^{8\pi i}/9$, $e^{14\pi i}/9$
// ZZ is the basic primitive 9'th root of unity. 
// ZZ2, ZZ4, and ZZ7 are powers of ZZ.
#define ZZ (0.766044443118978035202392650555+0.642787609686539326322643409907*I)
#define ZZ2 (0.173648177666930348851716626769+0.984807753012208059366743024590*I)
#define ZZ4 (-0.939692620785908384054109277325+0.342020143325668733044099614682*I)
#define ZZ7 (0.173648177666930348851716626769-0.984807753012208059366743024590*I)
// WW = ZZ+ZZ^-1 = $2\cos(2\pi/9)$
#define WW (1.532088886237956070404785301111)

#define MaxNoGens 4
#define ActualNoGens 4
#define KSiz 1

void SetupGens(Elt_t GenElts[],int *NoGenElts_p);
void SetupK(Elt_t KElts[]);

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 GenIx,GenEltIx,NoGenElts;
  size_t NamLen;
  Mtx_t MtxInv;

  FILE *GensFile;
  FILE *KFile;

  // 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);
  }

  // Prepare the generators and write the generator file
  SetupGens(GenElts,&NoGenElts);
  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);
  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) {
  // The form with respect to which the original, given matrices are
  // unitary:
  Mtx_t FF0=
    {{ (3-3*sq6) + sq6*WW + sq6*WW*WW, 0, 0},
     { 0, (3+sq6) - sq6*WW*WW, 0},
     { 0, 0, (3-sq6) - sq6*WW}};

  // Matrix to conjugate them to be unitary with respect to
  // $\diag{1,1,-1}$.  The nonzero elements here are
  // obtained from the diagonal elements of _FF0_ by taking square
  // roots of absolute values.
  Mtx_t ConjMtxInv,ConjMtx={
    {0,0,1.789504332114106653048452267768},
    {0,0.547894700539667349646553576350,0},
    {2.270252217086258544077485697624,0,0}};

  // The generators
  Gen_t Gens[MaxNoGens]={
    {"B","b",20,
     {
       {1/36.0*((sq6*sqm3-3*sq6-12)*ZZ2+(-4*sq6*sqm3-6*sq6-24)*ZZ),
	1/36.0*(((2*sq6+18)*sqm3-6)*ZZ2+(-sq6*sqm3+(3*sq6-24))*ZZ
		+2*sq6*sqm3+12*sq6+12),
	1/36.0*(((sq6+12)*sqm3+(3*sq6+12))*ZZ2+(-5*sq6*sqm3-3*sq6-24)*ZZ
		-2*sq6*sqm3+12*sq6+12)},
       {1/36.0*((-4*sq6*sqm3+(6*sq6+12))*ZZ2+((-sq6+6)*sqm3+(3*sq6+30))*ZZ
		+4*sq6*sqm3+6*sq6+24),
	1/36.0*(((sq6+6)*sqm3+(3*sq6+6))*ZZ2+((-sq6-12)*sqm3+(9*sq6+12))*ZZ),
	1/36.0*(((-sq6-6)*sqm3+(3*sq6+30))*ZZ2+((2*sq6-12)*sqm3+12)*ZZ
		+2*sq6*sqm3+12*sq6+12)},
       {1/36.0*(((-5*sq6-10)*sqm3-7*sq6+6)*ZZ2+((3*sq6+2)*sqm3+(5*sq6-18))*ZZ
		+(4*sq6+8)*sqm3+2*sq6+24),
	1/36.0*(((-sq6-6)*sqm3-9*sq6-6)*ZZ2+((2*sq6+12)*sqm3-24)*ZZ
		+4*sq6*sqm3+6*sq6+24),
	1/36.0*(((-2*sq6-6)*sqm3+6)*ZZ2+((5*sq6+12)*sqm3-3*sq6+12)*ZZ)}
     }
    },

    {"D","d",20,
     {
       {1/36.0*((3*sq6*sqm3+9*sq6)*ZZ2+(-sq6*sqm3+(3*sq6+12))*ZZ+(4*sq6+6)*sqm3+6),
	1/36.0*(((-7*sq6-12)*sqm3-3*sq6)*ZZ2+((5*sq6+18)*sqm3+(3*sq6+18))*ZZ
		+(-7*sq6-6)*sqm3-3*sq6-18),
	1/36.0*(((6*sq6+18)*sqm3-12*sq6-18)*ZZ2+((-8*sq6-12)*sqm3+(6*sq6+12))*ZZ
		+(5*sq6+24)*sqm3-3*sq6-12)},
       {1/36.0*(((5*sq6+18)*sqm3+(9*sq6+30))*ZZ2+((7*sq6+18)*sqm3-3*sq6+6)*ZZ
		+(7*sq6+6)*sqm3-9*sq6-30),
	1/36.0*(-6*sq6*sqm3*ZZ2+((2*sq6+6)*sqm3-6)*ZZ+(4*sq6+6)*sqm3+6),
	1/36.0*(((5*sq6+6)*sqm3-9*sq6-18)*ZZ2+(-sq6*sqm3-9*sq6-36)*ZZ
		+(-7*sq6-6)*sqm3-3*sq6-18)},
       {1/36.0*(((6*sq6+10)*sqm3+(4*sq6+18))*ZZ2+((-4*sq6-2)*sqm3+(10*sq6+30))*ZZ
		+(-3*sq6-20)*sqm3+sq6),
	1/36.0*(((-7*sq6-24)*sqm3+(3*sq6+12))*ZZ2+((-5*sq6-6)*sqm3-9*sq6-30)*ZZ
		+(7*sq6+6)*sqm3-9*sq6-30),
	1/36.0*((3*sq6*sqm3-9*sq6)*ZZ2+((-sq6-6)*sqm3-3*sq6-6)*ZZ+(4*sq6+6)*sqm3+6)}
     }
    },

    {"E","e",20,
     {
       {1/18.0*((sq6*sqm3-3*sq6-12)*ZZ2+((sq6+3)*sqm3-3*sq6-3)*ZZ
		+(sq6+3)*sqm3-3*sq6-3),
	1/18.0*(((-sq6-6)*sqm3+(3*sq6+6))*ZZ2+((2*sq6-3)*sqm3-6*sq6-3)*ZZ
		-sq6*sqm3-3*sq6+6),
	1/18.0*(((3*sq6+3)*sqm3+(3*sq6+9))*ZZ2+((-sq6-3)*sqm3-3*sq6-3)*ZZ
		+(sq6+6)*sqm3+3*sq6+12)},
       {1/18.0*(((-sq6-6)*sqm3+(3*sq6+12))*ZZ2+(sq6*sqm3+(3*sq6+6))*ZZ
		+(3*sq6+6)*sqm3+3*sq6),
	1/18.0*(((sq6+6)*sqm3+(3*sq6+6))*ZZ2+((-2*sq6-3)*sqm3-3)*ZZ
		+(sq6+3)*sqm3-3*sq6-3),
	1/18.0*((-sq6*sqm3-3*sq6-12)*ZZ2+(-4*sq6*sqm3+6)*ZZ-sq6*sqm3-3*sq6+6)},
       {1/18.0*(((2*sq6+6)*sqm3-6)*ZZ2+((3*sq6+3)*sqm3-3*sq6+9)*ZZ-sq6*sqm3+3*sq6-6),
	1/18.0*(((-sq6-3)*sqm3-3*sq6-15)*ZZ2+((sq6+3)*sqm3-3*sq6-3)*ZZ
		+(3*sq6+6)*sqm3+3*sq6),
	1/18.0*(((-2*sq6-6)*sqm3+6)*ZZ2+(sq6*sqm3+(3*sq6+6))*ZZ+(sq6+3)*sqm3-3*sq6-3)}
     }
    },

    {"U","u",20,
     {
       {(-1/36.0*sq6*sqm3+1/12.0*(sq6+4))*ZZ2+(1/18.0*(sq6+3)*sqm3-1/6.0)*ZZ
	+1/18.0*(sq6+3)*sqm3+1/6.0*(-sq6-1),
	(1/9.0*(-sq6-3)*sqm3+1/6.0*(sq6-2))*ZZ2
	+(1/36.0*(11*sq6+6)*sqm3+1/12.0*(sq6+2))*ZZ
	+1/18.0*(-2*sq6-9)*sqm3+1/6.0*(-sq6+1),
	(1/12.0*(sq6+6)*sqm3+1/12.0*(-sq6-6))*ZZ2
	+(1/36.0*(-7*sq6-12)*sqm3+1/12.0*(sq6-4))*ZZ
	+1/18.0*(2*sq6+3)*sqm3+1/6.0*(sq6-1)},
       {(1/18.0*(sq6+6)*sqm3+1/3.0*(sq6+1))*ZZ2
	+(1/12.0*(sq6+6)*sqm3+1/12.0*(sq6+6))*ZZ
	+1/18.0*(2*sq6+3)*sqm3+1/6.0*(-sq6+1),
	(1/36.0*(-sq6-6)*sqm3+1/12.0*(-sq6-2))*ZZ2
	+(1/36.0*(-sq6-6)*sqm3+1/12.0*(-sq6-2))*ZZ
	+1/18.0*(sq6+3)*sqm3+1/6.0*(-sq6-1),
	(1/36.0*(-sq6+12)*sqm3+1/12.0*(-3*sq6-4))*ZZ2
	+(-1/9.0*sq6*sqm3+1/6.0*(-3*sq6-2))*ZZ
	+1/18.0*(-2*sq6-9)*sqm3+1/6.0*(-sq6+1)},
       {(1/36.0*(sq6+18)*sqm3+1/12.0*(sq6+2))*ZZ2
	+(1/12.0*(sq6-4)*sqm3+1/4.0*(sq6+4))*ZZ
	+1/18.0*(-4*sq6-3)*sqm3+1/6.0*(sq6+1),
	(1/36.0*(-7*sq6-12)*sqm3+1/12.0*(-sq6+4))*ZZ2+1/6.0*(-sq6-6)*ZZ
	+1/18.0*(2*sq6+3)*sqm3+1/6.0*(-sq6+1),
	(1/18.0*(sq6+3)*sqm3-1/6.0)*ZZ2+(-1/36.0*sq6*sqm3+1/12.0*(sq6+4))*ZZ
	+1/18.0*(sq6+3)*sqm3+1/6.0*(-sq6-1)}
     }
    }
  };

  int GenIx,GenEltIx;
  size_t NamLen;
  Mtx_t MtxInv;
  
  // Calculate the inverse of ConjMtx
  InvMtx(ConjMtxInv,ConjMtx);
  
  GenEltIx=0;
  for(GenIx=0;GenIx<ActualNoGens;GenIx++) {
    printf("Checking Generator %s\n",Gens[GenIx].Nam);
    // Check that the matrix is unitary with respect to FF0
    assert(CheckUnitary(Gens[GenIx].Mtx,FF0));
    // Calculate the conjugated matrix
    ConjugateMtx(GenElts[GenEltIx].Mtx,Gens[GenIx].Mtx,ConjMtx,ConjMtxInv);
    // and check that that is unitary relative to FF1
    assert(CheckUnitary(GenElts[GenEltIx].Mtx,FF1));
    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(CheckUnitary(GenElts[GenEltIx].Mtx,FF1));
      // 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 (trivial) finite group, K
//
void SetupK(Elt_t KElts[]) {
  // The identity element
  Gen_t IGen=
    {"1","",0,{{1,0,0},{0,1,0},0,0,1}};
  Elt_t *KElt_p;
  
  KElt_p=&KElts[0];
  
  // Set up the matrix
  memcpy(KElt_p->Mtx,IGen.Mtx,sizeof(Mtx_t));
  // the generator name
  assert(strlen(IGen.Nam)<MaxGenNameLength);
  strcpy(KElt_p->Word,IGen.Nam);
  // the ``length''
  KElt_p->Len = IGen.Len;
  // and KIx and KStart
  KElt_p->KIx = 0;
  KElt_p->KStart = 0;
  
  FixUpElt(KElt_p);
}