#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 "C10p2/C10p2-17"

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{2}$ 
#define sq2 (1.414213562373095048801688724209)
// $U = ((1+\sqrt{2}) + \sqrt{-5+2\sqrt{2}}) / 2$
#define UU (1.207106781186547524400844362105+0.736812879103950296577504567300*I)
// WW = $2\cos(2\pi/9)$
#define WW (1.532088886237956070404785301111)
// WW2 = WW^2
#define WW2 (2.347296355333860697703433253539)

#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=
    {{ -sq2+(1-sq2)*WW+WW2, 0, 0},
     { 0, (4-3*sq2)+sq2*WW+(-1+sq2)*WW2, 0},
     { 0, 0, (2+sq2)-WW-sq2*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.546324809540912490675299134410,0,0},
    {0,1.701864324747346906543242431920,0},
    {0,0,1.198938557427753470429266518300}};

  // The generators
  Gen_t Gens[MaxNoGens]={
    {"A","a",30,
     {
       {1/153.0*(((24*sq2+26)*UU-59*sq2-71)*WW2+((-30*sq2-41)*UU+(44*sq2+59))*WW+(-99*sq2-103)*UU+169*sq2+193),
	1/153.0*(((26*sq2+31)*UU-37*sq2-33)*WW2+((58*sq2+77)*UU-59*sq2-105)*WW+(-31*sq2-35)*UU+50*sq2+57),
	1/153.0*(((-15*sq2-63)*UU+(73*sq2+21))*WW2+((-42*sq2-3)*UU-37*sq2+69)*WW+(60*sq2+99)*UU-122*sq2-237)},
       {1/153.0*(((9*sq2+82)*UU-3*sq2-84)*WW2+((-27*sq2+77)*UU+(60*sq2-54))*WW+(-99*sq2-137)*UU+33*sq2+108),
	1/153.0*(((30*sq2+41)*UU-44*sq2-59)*WW2+((54*sq2+67)*UU-103*sq2-130)*WW+(-111*sq2-133)*UU+139*sq2+169),
	1/153.0*(((-58*sq2-77)*UU+(59*sq2+105))*WW2+((-32*sq2-46)*UU+(22*sq2+72))*WW+(137*sq2+181)*UU-142*sq2-219)},
       {1/153.0*(((19*sq2+56)*UU-46*sq2-64)*WW2+((-25*sq2-20)*UU+(31*sq2+52))*WW+(-8*sq2-88)*UU+65*sq2+86),
	1/153.0*(((27*sq2-77)*UU-60*sq2+54)*WW2+((36*sq2+5)*UU-63*sq2-30)*WW+(-135*sq2+181)*UU+147*sq2-168),
	1/153.0*(((-54*sq2-67)*UU+(103*sq2+130))*WW2+((-24*sq2-26)*UU+(59*sq2+71))*WW+(57*sq2+83)*UU-155*sq2-209)}
     }
    },

    {"B","b",20,
     {
       {1/306.0*(((-42*sq2-122)*UU+(48*sq2+52))*WW2+((-24*sq2-94)*UU+(42*sq2+20))*WW+(186*sq2+346)*UU-198*sq2-308),
	1/306.0*(((82*sq2+120)*UU-118*sq2-6)*WW2+((128*sq2+252)*UU-122*sq2-288)*WW+(16*sq2-96)*UU+74*sq2-138),
	1/306.0*(((22*sq2+4)*UU+(38*sq2-24))*WW2+((14*sq2+86)*UU-50*sq2-6)*WW+(82*sq2+52)*UU-118*sq2-108)},
       {1/306.0*(((-54*sq2-50)*UU+(86*sq2+28))*WW2+((-42*sq2+14)*UU+(82*sq2-16))*WW+(186*sq2+244)*UU-232*sq2-308),
	1/306.0*(((24*sq2+94)*UU-42*sq2-20)*WW2+((-18*sq2-28)*UU+(6*sq2+32))*WW+(54*sq2-86)*UU-18*sq2-164),
	1/306.0*(((-128*sq2-252)*UU+(122*sq2+288))*WW2+((-46*sq2-132)*UU+(4*sq2+282))*WW+(436*sq2+648)*UU-406*sq2-726)},
       {1/306.0*(((-29*sq2+174)*UU-132*sq2-92)*WW2+((-139*sq2-84)*UU+(120*sq2+164))*WW+(235*sq2-186)*UU+120*sq2-176),
	1/306.0*(((42*sq2-14)*UU-82*sq2+16)*WW2+((-12*sq2-64)*UU+(4*sq2+44))*WW+(-6*sq2+172)*UU+104*sq2-284),
	1/306.0*(((18*sq2+28)*UU-6*sq2-32)*WW2+((42*sq2+122)*UU-48*sq2-52)*WW+(66*sq2+46)*UU-90*sq2-140)}
     }
    },

    {"C","c",50,
     {
       {1/153.0*(((5*sq2+4)*UU+(21*sq2+44))*WW2+((-2*sq2-22)*UU+(12*sq2-38))*WW+(41*sq2+94)*UU-144*sq2-241),
	1/153.0*(((72*sq2+61)*UU-58*sq2-111)*WW2+((90*sq2+89)*UU-98*sq2-177)*WW+(-78*sq2+43)*UU-25*sq2+99),
	1/153.0*(((-sq2+6)*UU+(23*sq2+32))*WW2+((55*sq2+78)*UU-41*sq2-128)*WW+(89*sq2+129)*UU-211*sq2-247)},
       {1/153.0*(((-30*sq2-92)*UU+(78*sq2+110))*WW2+((-12*sq2-64)*UU+(72*sq2+112))*WW+(126*sq2+247)*UU-297*sq2-394),
	1/153.0*(((2*sq2+22)*UU-12*sq2+38)*WW2+((7*sq2+26)*UU+(9*sq2+82))*WW+(47*sq2+58)*UU-78*sq2-229),
	1/153.0*(((-90*sq2-89)*UU+(98*sq2+177))*WW2+((-18*sq2-28)*UU+(40*sq2+66))*WW+(246*sq2+343)*UU-337*sq2-477)},
       {1/153.0*(((-sq2+6)*UU-28*sq2-36)*WW2+((-47*sq2-75)*UU+(61*sq2+144))*WW+(89*sq2+78)*UU-109*sq2-60),
	1/153.0*(((12*sq2+64)*UU-72*sq2-112)*WW2+((-18*sq2-28)*UU+(6*sq2-2))*WW+(42*sq2-65)*UU+3*sq2+50),
	1/153.0*(((-7*sq2-26)*UU-9*sq2-82)*WW2+((-5*sq2-4)*UU-21*sq2-44)*WW+(65*sq2+154)*UU-84*sq2+11)}
     }
    }
  };

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