#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 "a1p5/a1p5-0"

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{5}$ 
#define sq5 (2.236067977499789696409173668731)
// Root of Z^3-3Z^2-2
#define ZZ (3.195823345445647152832799205550)
// $\sqrt{2}$
#define sq2 (1.414213562373095048801688724210)

#define MaxNoGens 3
#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=
    {{ 5, 0, 0},
     { 0, 0, 1},
     { 0, 1, 0}};

  // Matrix to conjugate them to be unitary with respect to
  // $\diag{1,1,-1}$
  Mtx_t ConjMtx={{sq5,0,0},{0,1/sq2,1/sq2},{0,-1/sq2,1/sq2}},ConjMtxInv;

  // The generators
  Gen_t Gens[MaxNoGens]={
    {"B","b",2,
     {
       {1/18.0*(-2*I+1)*ZZ*ZZ+1/18.0*(11*I+2)*ZZ+1/18.0*(-5*I-5),
	1/18.0*(-I+2)*ZZ*ZZ+1/18.0*(-2*I-5)*ZZ+1/18.0*(11*I-1),
	1/90.0*(-7*I-4)*ZZ*ZZ+1/90.0*(16*I+7)*ZZ+1/90.0*(11*I+17)
       },
       {1/18.0*(-I+8)*ZZ*ZZ+1/18.0*(4*I-17)*ZZ+1/18.0*(5*I+5),
	1/9.0*(2*I-1)*ZZ*ZZ+1/9.0*(-8*I+4)*ZZ+1/9.0*(2*I-1),
	1/9.0*(I+1)*ZZ*ZZ+1/9.0*(-I-4)*ZZ+1/9.0*(I+1)
       },
       {1/18.0*(-5*I-20)*ZZ*ZZ+1/18.0*(20*I+65)*ZZ+1/18.0*(25*I-5),
	1/9.0*(4*I-2)*ZZ*ZZ+1/9.0*(-10*I+5)*ZZ+1/9.0*(I+7),
	1/18.0*(-2*I+1)*ZZ*ZZ+1/18.0*(5*I-10)*ZZ+1/18.0*(I+7)
       }
     },
    },
    {"A","a",10,
      {
	{1/18.0*(-7*I+3)*ZZ*ZZ+1/18.0*(19*I+3)*ZZ+1/9.0*(-2*I-12),
	 1/18.0*I*ZZ*ZZ+1/18.0*(-4*I-3)*ZZ+1/18.0*(I+9),
	 1/9.0*I*ZZ*ZZ-1/9.0*I*ZZ-2/9.0*I
	},
	{1/18.0*(-5*I-15)*ZZ*ZZ+1/18.0*(5*I+45)*ZZ+5/9.0*I,
	 1/18.0*(11*I+3)*ZZ*ZZ+1/18.0*(-35*I-9)*ZZ+1/9.0*(-2*I-6),
	 -1/9.0*I*ZZ*ZZ+4/9.0*I*ZZ+1/9.0*(-I+3)
	},
	{5/18.0*I*ZZ*ZZ+1/18.0*(-20*I+15)*ZZ+1/18.0*(5*I+15),
	 1/18.0*(-5*I+15)*ZZ*ZZ+1/18.0*(5*I-45)*ZZ+5/9.0*I,
	 1/9.0*(-2*I-3)*ZZ*ZZ+1/9.0*(8*I+3)*ZZ-5/9.0*I
	}
      },
    },
    {"C","c",15,
     {
       {1/18.0*(-3*I-2)*ZZ*ZZ-1/18.0*ZZ+1/18.0*(3*I+13),
	1/18.0*(6*I+5)*ZZ*ZZ+1/18.0*(-21*I-14)*ZZ+1/18.0*(15*I-7),
	1/18.0*ZZ*ZZ+1/18.0*(-3*I-4)*ZZ+1/18.0*(3*I-5)
       },
       {5/18.0*ZZ*ZZ+1/18.0*(15*I-20)*ZZ+1/18.0*(-15*I-25),
	1/18.0*(-3*I+7)*ZZ*ZZ+1/18.0*(15*I-19)*ZZ+1/9.0*(-6*I+2),
	-2/9.0*ZZ*ZZ+1/9.0*(3*I+2)*ZZ+1/9.0*(3*I+1)
       },
       1/18.0*(-30*I-5)*ZZ*ZZ+1/18.0*(75*I+50)*ZZ+1/18.0*(75*I-65),
       -5/9.0*ZZ*ZZ+20/9.0*ZZ-20/9.0,
       1/18.0*(6*I-5)*ZZ*ZZ+1/18.0*(-15*I+20)*ZZ+1/18.0*(-9*I+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++) {
    // 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);
}