#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 "C20p2/C20p2-3m"

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{7}$ 
#define sq7 (2.6457513110645905905016158)
// $\exp(2\pi i/7)$, a seventh root of unity
#define ZZ (0.6234898018587335305300629+0.7818314824680298087044108*I)
// and its powers, the other seventh roots of unity
#define ZZ2 (ZZ*ZZ)
#define ZZ3 (ZZ*ZZ2)
#define ZZ4 (ZZ*ZZ3)
#define ZZ5 (ZZ*ZZ4)
#define ZZ6 (ZZ*ZZ5)
// WW = ZZ+ZZ^-1
#define WW (1.2469796037174670610500098)
#define WW2 (WW*WW)

#define MaxNoGens 5
#define ActualNoGens 3
#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=
    {{(2-4*sq7/7.0) + (1-3*sq7/7.0)*WW + (-sq7/7.0)*WW2, 0, 0},
     { 0, (-sq7/7.0) + (sq7/7.0)*WW + (1-2*sq7/7.0)*WW2, 0},
     { 0, 0, (3-9*sq7/7.0) + (-1+2*sq7/7.0)*WW + (-1+3*sq7/7.0)*WW2}};

  // 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.5162743317096442384422437,0,0},
    {0,0,0.7055729746140944909267776},
    {0,0.6876552813429416382391305,0}};

  Gen_t Gens[MaxNoGens]={
    {"B","b",10,
     {
       {(1/28.0)*(((-2*sq7-14)*I-6*sq7-14)*ZZ2+(-2*sq7*I-8*sq7-14)*ZZ+(4*sq7+28)*I),
	(1/28.0)*(-4*sq7*ZZ2+(-4*sq7*I+4*sq7)*ZZ+4*sq7*I),
	(1/28.0)*(((-6*sq7-14)*I+(2*sq7+14))*ZZ2+(-14*I-2*sq7)*ZZ-8*sq7*I)},
       {(1/28.0)*(((sq7-7)*I-3*sq7-21)*ZZ2+((sq7+7)*I-5*sq7-21)*ZZ+(5*sq7+21)*I+sq7+7),
	(1/28.0)*((14*I-2*sq7)*ZZ2+((2*sq7+14)*I+(6*sq7+14))*ZZ+(-2*sq7+14)*I+10*sq7+14),
	(1/28.0)*((-4*sq7*I+8*sq7)*ZZ2+4*sq7*ZZ+(4*sq7-14)*I+2*sq7)},
       {(1/28.0)*(((3*sq7-7)*I-3*sq7-7)*ZZ2+((3*sq7-14)*I-6*sq7-7)*ZZ+(8*sq7+7)*I-5*sq7),
	(1/28.0)*((14*I-2*sq7)*ZZ2+((-sq7+7)*I+(3*sq7+21))*ZZ+(-6*sq7+14)*I+6*sq7+14),
	(1/28.0)*((2*sq7*I+(8*sq7+14))*ZZ2+(-14*I+2*sq7)*ZZ-2*sq7*I+4*sq7+14)}
     }
    },

    {"A","a",10,
     {
       {(1/84.0)*((4*sq7+56)*I*ZZ2+((8*sq7+28)*I+12*sq7)*ZZ+(2*sq7-14)*I+30*sq7+42),
	(1/84.0)*(((-4*sq7-14)*I-10*sq7+28)*ZZ2+(-28*sq7+28)*I*ZZ+(4*sq7+14)*I+10*sq7-28),
	(1/84.0)*(((4*sq7+56)*I+(8*sq7+28))*ZZ2+((3*sq7+21)*I+(15*sq7+63))*ZZ+(-14*sq7+14)*I+26*sq7-14)},
       {(1/84.0)*(((8*sq7-14)*I+(14*sq7+28))*ZZ2+((-10*sq7-56)*I+(8*sq7-14))*ZZ-12*sq7*I-8*sq7-28),
	(1/84.0)*(((4*sq7-28)*I+12*sq7)*ZZ2+(-4*sq7-56)*I*ZZ-42*I+2*sq7+28),
	(1/84.0)*(((-24*sq7+42)*I+(10*sq7-28))*ZZ2+((4*sq7+14)*I+(10*sq7-28))*ZZ+(20*sq7-14)*I+22*sq7-28)},
       {(1/84.0)*(((21*sq7-21)*I+(7*sq7-49))*ZZ2+((25*sq7-49)*I+(3*sq7-21))*ZZ+(31*sq7-7)*I-3*sq7+21),
	(1/84.0)*(((-18*sq7-42)*I-6*sq7-42)*ZZ2+((-8*sq7+14)*I-14*sq7-28)*ZZ+(-2*sq7+56)*I-8*sq7-70),
	(1/84.0)*(((-8*sq7-28)*I-12*sq7)*ZZ2+((-4*sq7+28)*I-12*sq7)*ZZ+(-2*sq7+14)*I+10*sq7+14)}
     }
    },

    {"C","c",10,
     {
       {(1/84.0)*(((-8*sq7-28)*I-16*sq7-56)*ZZ2+((8*sq7-14)*I-2*sq7-28)*ZZ+(14*sq7+28)*I+4*sq7+14),
	(1/84.0)*(((-20*sq7-28)*I-12*sq7+84)*ZZ2+((-40*sq7+28)*I-8*sq7-28)*ZZ+(-10*sq7-14)*I+6*sq7-42),
	(1/84.0)*(((6*sq7+42)*I-2*sq7+14)*ZZ2+(42*I+6*sq7)*ZZ+(15*sq7+21)*I+31*sq7+35)},
       {(1/84.0)*(((-8*sq7+14)*I+6*sq7)*ZZ2+((-8*sq7-28)*I+12*sq7)*ZZ+(16*sq7+14)*I+10*sq7-28),
	(1/84.0)*(((16*sq7+14)*I+(14*sq7+28))*ZZ2+((8*sq7+28)*I+(16*sq7+56))*ZZ+(-10*sq7-14)*I+26*sq7+70),
	(1/84.0)*(((-20*sq7+56)*I+(4*sq7-112))*ZZ2+((20*sq7+28)*I+(12*sq7-84))*ZZ+(42*sq7-42)*I+26*sq7-14)},
       {(1/84.0)*(((37*sq7-7)*I+(13*sq7-49))*ZZ2+((43*sq7-49)*I+(11*sq7+49))*ZZ+(25*sq7-49)*I+11*sq7+91),
	(1/84.0)*((-42*I+6*sq7)*ZZ2+((8*sq7-14)*I-6*sq7)*ZZ+(20*sq7+28)*I),
	(1/84.0)*(((-8*sq7+14)*I+(2*sq7+28))*ZZ2+((-16*sq7-14)*I-14*sq7-28)*ZZ+(-4*sq7+28)*I+12*sq7)}
     }
    },

    {"A20","a20",10,
     {
       {(1/84.0)*(((8*sq7+28)*I-40*sq7-56)*ZZ2+((-8*sq7+56)*I-20*sq7-28)*ZZ+(28*sq7+140)*I+4*sq7-28),
	(1/84.0)*((12*sq7*I+(8*sq7+28))*ZZ2+((-14*sq7-28)*I-8*sq7+14)*ZZ+(-12*sq7-42)*I-14*sq7-28),
	(1/84.0)*(((34*sq7+14)*I-6*sq7+42)*ZZ2+((-3*sq7+21)*I+(15*sq7+105))*ZZ+(-38*sq7-28)*I+12*sq7+42)},
       {(1/84.0)*(((-12*sq7+42)*I+(14*sq7+112))*ZZ2+(-36*sq7*I+(8*sq7+28))*ZZ+(-50*sq7-28)*I+20*sq7-14),
	(1/84.0)*(((-16*sq7+28)*I+(20*sq7+28))*ZZ2+((-8*sq7-28)*I+(40*sq7+56))*ZZ+(-4*sq7-14)*I+10*sq7-28),
	(1/84.0)*(((-26*sq7-28)*I-16*sq7-14)*ZZ2+(-12*sq7*I-8*sq7-28)*ZZ+(-4*sq7-14)*I-18*sq7-84)},
       {(1/84.0)*(((7*sq7+35)*I-sq7-35)*ZZ2+((23*sq7+49)*I+(5*sq7-35))*ZZ+(5*sq7-35)*I+17*sq7+7),
	(1/84.0)*(((-24*sq7-42)*I-6*sq7-84)*ZZ2+((12*sq7-42)*I-14*sq7-112)*ZZ+12*sq7*I-8*sq7-28),
	(1/84.0)*(((8*sq7-56)*I+(20*sq7+28))*ZZ2+((16*sq7-28)*I-20*sq7-28)*ZZ+(18*sq7+42)*I-14*sq7+14)}
     }
    },

    {"B2","b2",10,
     {
       {(1/28.0)*((-2*sq7*I+(4*sq7+14))*ZZ2+2*sq7*I-4*sq7-14),
	(1/28.0)*(((6*sq7+28)*I+14)*ZZ2+((4*sq7+14)*I+10*sq7)*ZZ+4*sq7*I+4*sq7+28),
	(1/28.0)*(((-8*sq7-14)*I-2*sq7)*ZZ2+(-6*sq7*I-8*sq7-14)*ZZ+-4*sq7-28)},
       {(1/28.0)*((-2*sq7*I-8*sq7-14)*ZZ2+((6*sq7+14)*I-2*sq7-14)*ZZ+(3*sq7+7)*I+3*sq7-7),
	(1/28.0)*((2*sq7*I-4*sq7-14)*ZZ2+(2*sq7*I-4*sq7-14)*ZZ+(10*sq7+14)*I-6*sq7-14),
	(1/28.0)*(((-2*sq7-14)*I+(10*sq7-14))*ZZ2+((-6*sq7-28)*I-14)*ZZ+(8*sq7-14)*I-10*sq7)},
       {(1/28.0)*(((-3*sq7-28)*I-4*sq7+7)*ZZ2+((5*sq7-7)*I-11*sq7+7)*ZZ+(5*sq7+7)*I-13*sq7+7),
	(1/28.0)*(((8*sq7+14)*I+6*sq7)*ZZ2+(2*sq7*I+(8*sq7+14))*ZZ+(-3*sq7-21)*I+7*sq7+7),
	(1/28.0)*((-2*sq7*I+(4*sq7+14))*ZZ+2*sq7*I-4*sq7-14)}
     }
    }
  };

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