// THIS IS UNFINISHED

#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-0-KK"

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)

#define MaxNoGens 2
#define ActualNoGens 2
#define KSiz 7

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 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,FF0,ConjMtx,ConjMtxInv);
  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,ConjMtx,ConjMtxInv);
  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 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]={
    {"A","a",20,
     {
       {(1/7.0*(-3*sq7-7)*I-1/7.0*sq7)*ZZ2+(1/14.0*(-3*sq7-7)*I+1/14.0*(-sq7-7))*ZZ+1/7.0*sq7*I+1/7.0*(-2*sq7-7),
	(1/14.0*(sq7+7)*I+1/14.0*(3*sq7+7))*ZZ2+(1/14.0*(-5*sq7-7)*I+1/2.0*(sq7+1))*ZZ+1/14.0*(-3*sq7-14)*I+1/14.0*(4*sq7+7),
	(1/14.0*(-2*sq7+7)*I+1/2.0*(-sq7-2))*ZZ2+(1/14.0*(2*sq7+7)*I+1/14.0*(-3*sq7-14))*ZZ+1/4.0*(sq7+3)*I+1/28.0*(-sq7-21)},
       {(1/14.0*(3*sq7+7)*I+1/14.0*(3*sq7-7))*ZZ2+(1/14.0*(5*sq7+7)*I+1/2.0*(sq7+1))*ZZ+1/14.0*(-sq7-14)*I+1/14.0*(4*sq7+7),
	(1/14.0*(3*sq7+7)*I+1/14.0*(sq7-7))*ZZ2+(1/7.0*(3*sq7+7)*I+1/7.0*sq7)*ZZ+5/14.0*sq7*I+1/14.0*(4*sq7+7),
	(1/7.0*(-3*sq7-7)*I+2/7.0*sq7)*ZZ2+(1/14.0*(-sq7-7)*I+1/14.0*(-3*sq7-7))*ZZ+-1/2.0*I-1/14.0*sq7},
       {(1/14.0*(2*sq7-7)*I+1/2.0*(-sq7-2))*ZZ2+(1/28.0*(9*sq7+7)*I+1/28.0*(-13*sq7-21))*ZZ+1/28.0*(15*sq7+35)*I+1/28.0*(-sq7+21),
	(1/7.0*sq7*I+1/7.0*(2*sq7+7))*ZZ2+(1/14.0*(-3*sq7-7)*I+1/14.0*(-3*sq7+7))*ZZ+1/14.0*(-6*sq7-7)*I-1/14.0*sq7,
	(1/14.0*(3*sq7+7)*I+1/14.0*(sq7+7))*ZZ2+(1/14.0*(-3*sq7-7)*I+1/14.0*(-sq7+7))*ZZ}
     }
    },

    {"C","c",20,
     {
       {(1/2.0*I+1/14.0*(5*sq7+14))*ZZ2+(1/14.0*(-3*sq7-7)*I+1/14.0*(sq7+7))*ZZ+1/14.0*(-4*sq7-7)*I+1/14.0*sq7,
	-1/7.0*sq7*ZZ2+(-1/7.0*sq7*I+1/7.0*sq7)*ZZ+1/7.0*sq7*I,
	(1/2.0*I-3/14.0*sq7)*ZZ2+(1/14.0*(sq7+14)*I+1/2.0)*ZZ+1/14.0*(-sq7+7)*I+1/14.0*(3*sq7-7)},
       {(1/28.0*(5*sq7+21)*I+1/28.0*(sq7+7))*ZZ2+2/7.0*sq7*ZZ+1/28.0*(-5*sq7-7)*I+1/28.0*(5*sq7-7),
	(1/14.0*(-3*sq7-14)*I+1/14.0*(-4*sq7-7))*ZZ2+(-1/2.0*I+1/14.0*(-5*sq7-14))*ZZ+1/14.0*(3*sq7+7)*I+1/14.0*(-5*sq7-7),
	(-1/7.0*sq7*I+2/7.0*sq7)*ZZ2+1/7.0*sq7*ZZ+1/14.0*(2*sq7-7)*I+1/14.0*sq7},
       {(1/28.0*(3*sq7-7)*I+1/28.0*(-3*sq7-7))*ZZ2+(1/28.0*(3*sq7-14)*I+1/28.0*(-6*sq7-7))*ZZ+1/28.0*(8*sq7+7)*I-5/28.0*sq7,
	(1/28.0*(-5*sq7-21)*I+1/4.0*(sq7-1))*ZZ2+(1/28.0*(-5*sq7-21)*I+1/28.0*(-sq7-7))*ZZ+1/14.0*(-2*sq7-7)*I+1/14.0*(-3*sq7-14),
	(1/14.0*(3*sq7+7)*I+1/14.0*(-sq7-7))*ZZ2+(1/14.0*(3*sq7+14)*I+1/14.0*(4*sq7+7))*ZZ+1/14.0*sq7*I+1/14.0*(4*sq7+7)}
     }
    }
  };

  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 7-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 KMtx;
  // Powers of
  //
  //   ZZ = \zeta_7 = e^{2\pi i/7}
  //           
  COMPLEX_t Zeta7Pow[7]={1,ZZ,ZZ2,ZZ3,ZZ4,ZZ5,ZZ6};

  KElt_p=KElts;

  // Set up elements of the form:
  //
  //   ZZ^j     0        0
  //     0    ZZ^{2j}    0
  //     0      0      ZZ^{4j}
  //
  //   = Z^j

  
  KElt_p=&KElts[0];
  
  for(Ix1=0;Ix1<7;Ix1++) {
    // Set up the matrix
    Ix2=(2*Ix1) % 7;
    Ix3=(4*Ix1) % 7;
    memset(KElt_p->Mtx,0,sizeof(Mtx_t));
    KElt_p->Mtx[0][0]=Zeta7Pow[Ix1];
    KElt_p->Mtx[1][1]=Zeta7Pow[Ix2];
    KElt_p->Mtx[2][2]=Zeta7Pow[Ix3];
    if(Ix1==0) {
      // set up the generator name
      strncpy(KElt_p->Word,"1");
      // and the ``length''
      KElt_p->Len = 0;
    }
    else {
      // set up the generator name
      sprintf(KElt_p->Word,"Z^%i",Ix1);
      // and the ``length''
      KElt_p->Len = 2;
    }
    // Set up KIx and KStart
    KElt_p->KIx =  KElt_p->KStart = (KElt_p-KElts);

    FixUpElt(KElt_p);
    KElt_p++;
  }
}