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

// R0 =$\sqrt(\sqrt{3}-1))$
#define R0 (0.855599677167352192969235766210)

// The original generators act unitarily with respect to the
// form $\diag(1,1,1-\sqrt(3))=\diag(1,1,-R_0^2)$.

// Powers of $Z = -\zeta_{12}=e^{7\pi i/6}$
#define Z  (-0.866025403784438646763723170753-0.5*I)
#define Z2 ( 0.5+0.866025403784438646763723170753*I)
#define Z3 ( 0.0-1.0*I)
#define Z4 (-0.5+0.866025403784438646763723170753*I)
#define Z5 ( 0.866025403784438646763723170753-0.5*I)

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;

#define MaxNoGens 1
#define ActualNoGens 1
#define KSiz 288

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

  // The form with respect to which the original, given matrices are
  // unitary:
  Mtx_t FF0=
    {{ 1, 0,    0  },
     { 0, 1,    0  },
     { 0, 0, -R0*R0}};

  // Matrix by which to conjugate them to make them unitary with
  // respect to $\diag{1,1,-1}$
  Mtx_t ConjMtxInv, ConjMtx={
    { 1, 0,  0},
    { 0, 1,  0},
    { 0, 0, R0}};

  FILE *GensFile;
  FILE *KFile;

  // Direction in which to shift the origin:
  Pt_t Dir={0.35+0.45*I,0.55+0.65*I};
  // FLOAT_t sShift=0.600;   // How far to shift it (Euclidean distance)
  FLOAT_t sShift=0.550;   // How far to shift it (Euclidean distance)
  Pt_t W;

  Mtx_t C,CInv,ConjMtx2,ConjMtx2Inv;

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

  // Calculate the inverse of ConjMtx
  InvMtx(ConjMtxInv,ConjMtx);
  

  // Calculations for shifting the origin
  NormalizePt(Dir);
  for(Ix=0;Ix<2;Ix++) W[Ix]=sShift*Dir[Ix];
  OriginToZMtx(CInv,C,W); // C shifts W to 0; CInv shifts 0 to W
  
  MultMtx(ConjMtx2,C,ConjMtx);
  MultMtx(ConjMtx2Inv,ConjMtxInv,CInv);

  // Prepare the generators and write the generator file
  SetupGens(GenElts,&NoGenElts,FF0,ConjMtx2,ConjMtx2Inv);
  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,ConjMtx2,ConjMtx2Inv);
  fprintf(KFile,"KSiz=%i\n",KSiz,FF0,ConjMtx2,ConjMtx2Inv);
  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 generators
  Gen_t Gens[MaxNoGens]={
    {"A","a",100,
     {
       {1,0,0},
       {0,Z3-Z2,Z3-1},
       {0,-Z3+Z2+Z-1,-Z3+Z2-1}
     }
    },
  };

  int GenIx,GenEltIx;
  size_t NamLen;
  Mtx_t MtxInv,tmpMtx;

  GenEltIx=0;
  for(GenIx=0;GenIx<ActualNoGens;GenIx++) {
    // Check that the matrix is unitary with respect to FF0
    printf("Checking GenIx=%i\n",GenIx);
    // OutputMtx(stdout,Gens[GenIx].Mtx);
    assert(CheckPU(Gens[GenIx].Mtx,FF0,epsBig));
    // Calculate the conjugated matrix
    ConjugateMtx(GenElts[GenEltIx].Mtx,
		 Gens[GenIx].Mtx,ConjMtx,ConjMtxInv);
    // and check that that is unitary relative to FF1
    assert(CheckPU(GenElts[GenEltIx].Mtx,FF1,epsBigger));
    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(CheckPU(GenElts[GenEltIx].Mtx,FF1,epsBigger));
      // 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 288-element finite group, K
//
void SetupK(Elt_t KElts[],
	       Mtx_t FF0,Mtx_t ConjMtx,Mtx_t ConjMtxInv) {

  int Ix1,Ix2,RIx,CIx;
  Elt_t *KElt_p;
  Mtx_t tmpMtx,KMtx;
  // Powers of
  //
  //   -Z = \zeta_{12} = e^{\pi i/6}
  //           
  COMPLEX_t Zeta12Pow[12]
    ={1,-Z,Z2,-Z3,Z4,-Z5,-1,Z,-Z2,Z3,-Z4,Z5};

  // Exponent strings
  char ExpString[12][4]=
    {"??","","^2","^3","^4","^5","^6","^7","^8","^9","^10","^11"};
  char Word1[6],Word2[10];

  KElt_p=KElts;

  // Set up elements of the form:
  //
  //   (-Z)^{Ix1}    0        0
  //          0  (-Z)^{Ix2}   0
  //          0      0        1
  //
  //   = K1^{Ix2} K2 K1^{Ix1} K2

  memset(KMtx,0,sizeof(Mtx_t));
  KMtx[2][2]=1;

  for(Ix1=0;Ix1<12;Ix1++) {
    KMtx[0][0]=Zeta12Pow[Ix1];
    if(Ix1 == 0)
      strncpy(Word2,"",1);
    else {
      strncpy(&Word2[0],"K2K1",5);
      strncpy(&Word2[4],ExpString[Ix1],4);
      strncpy(&Word2[4+strlen(ExpString[Ix1])],"K2",3);
    }
    for(Ix2=0;Ix2<12;Ix2++) {
      KMtx[1][1]=Zeta12Pow[Ix2];
      if(Ix2 == 0) {
	if(Ix1 == 0)
	  strncpy(Word1,"1",2);
	else
	  strncpy(Word1,"",1);
      }
      else {
	strncpy(&Word1[0],"K1",3);
	strncpy(&Word1[2],ExpString[Ix2],4);
      }
      // Set up the matrix
      assert(CheckUnitary(KMtx,FF0));
      ConjugateMtx(KElt_p->Mtx,KMtx,ConjMtx,ConjMtxInv);
      assert(CheckPU(KElt_p->Mtx,FF1,epsSmall));
      // the generator name
      sprintf(KElt_p->Word,"%s%s",Word1,Word2);
      // the ``length''
      if(Ix1 == 0 && Ix2 ==0)
	KElt_p->Len = 0;
      else
	KElt_p->Len = 2;
      // KIx and KStart
      KElt_p->KIx = KElt_p->KStart = (KElt_p-KElts);

      FixUpElt(KElt_p);
      KElt_p++;
    }
  }

  // Set up elements of the form:
  //
  //          0    (-Z)^{Ix1}   0
  //   (-Z)^{Ix2}       0       0
  //          0         0        1
  //
  //   = K1^{Ix2} K2 K1^{Ix1}


  memset(KMtx,0,sizeof(Mtx_t));
  KMtx[2][2]=1;

  for(Ix1=0;Ix1<12;Ix1++) {
    KMtx[0][1]=Zeta12Pow[Ix1];
    if(Ix1 == 0)
      strncpy(Word2,"K2",3);
    else {
      strncpy(&Word2[0],"K2K1",5);
      strncpy(&Word2[4],ExpString[Ix1],4);
    }
    for(Ix2=0;Ix2<12;Ix2++) {
      KMtx[1][0]=Zeta12Pow[Ix2];
      if(Ix2 == 0) {
	strncpy(Word1,"",1);
      }
      else {
	strncpy(&Word1[0],"K1",3);
	strncpy(&Word1[2],ExpString[Ix2],4);
      }
      // Set up the matrix
      assert(CheckUnitary(KMtx,FF0));
      ConjugateMtx(KElt_p->Mtx,KMtx,ConjMtx,ConjMtxInv);
      assert(CheckPU(KElt_p->Mtx,FF1,epsSmall));
      // the generator name
      sprintf(KElt_p->Word,"%s%s",Word1,Word2);
      // the ``length''
      KElt_p->Len = 2;
      // KIx and KStart
      KElt_p->KIx = KElt_p->KStart = (KElt_p-KElts);

      FixUpElt(KElt_p);
      KElt_p++;
    }
  }
}