#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <complex.h>
#include <string.h>

FILE *EltsFile;
FILE *hsnorms;

#include "a7p2-0-base.c"

// Other generators of our group $\Gamma$
#define MaxNoGens00 3
#define MaxNoGens (2*MaxNoGens00)
#define NoGens00 2
#define NoGens (2*NoGens00)

Mtx_t Gens00[MaxNoGens00]={ 
  // These are the given generators
  // This is U'
  {
    {1.0/18*(-sqm7+3)*ZZ*ZZ+1.0/6*(-sqm7+1)*ZZ+1.0/6*(-sqm7-3),
     1.0/18*(sqm7+3)*ZZ*ZZ+1.0/6*(sqm7+3)*ZZ+1.0/3,
     1.0/126*(-sqm7-15)*ZZ*ZZ+1.0/42*(sqm7-19)*ZZ+1.0/21*(2*sqm7-2)},
    {1.0/63*(-sqm7+15)*ZZ*ZZ+1.0/21*(sqm7+19)*ZZ+1.0/21*(4*sqm7+4),
     1.0/9*sqm7*ZZ*ZZ+1.0/3*sqm7*ZZ+1.0/6*(-sqm7-1),
     1.0/126*(sqm7-21)*ZZ*ZZ+1.0/21*(-sqm7-7)*ZZ+1.0/14*(-sqm7+7)},
    {-8.0/63*sqm7*ZZ*ZZ+1.0/21*(-5*sqm7-7)*ZZ+1.0/21*(3*sqm7+7),
     2.0/63*sqm7*ZZ*ZZ-2.0/21*sqm7*ZZ-2.0/21*sqm7,
     1.0/18*(-sqm7-3)*ZZ*ZZ+1.0/6*(-sqm7-1)*ZZ+1.0/6*(-sqm7+1)}
  },

  // This is B'
  {
    {-2.0/63*sqm7*ZZ*ZZ+1.0/42*(sqm7+7)*ZZ+2.0/7*sqm7,
     -1.0/63*sqm7*ZZ*ZZ+1.0/42*(-3*sqm7-7)*ZZ+1.0/42*(-sqm7-7),
     1.0/126*(-sqm7+15)*ZZ*ZZ+1.0/84*(-sqm7+17)*ZZ+1.0/21*(sqm7-5)},
    {1.0/126*(13*sqm7-9)*ZZ*ZZ+1.0/42*(11*sqm7-3)*ZZ+1.0/42*(sqm7-1),
     1.0/126*(-sqm7+21)*ZZ*ZZ+1.0/42*(-5*sqm7+21)*ZZ+1.0/42*(3*sqm7-7),
     2.0/63*sqm7*ZZ*ZZ+1.0/7*sqm7*ZZ+1.0/21*sqm7},
    {-2.0/63*sqm7*ZZ*ZZ+1.0/21*(-3*sqm7+7)*ZZ+1.0/21*(-sqm7+7),
     1.0/126*(-11*sqm7-21)*ZZ*ZZ+1.0/21*(-5*sqm7-7)*ZZ+2.0/21*sqm7,
     1.0/126*(5*sqm7-21)*ZZ*ZZ+1.0/21*(2*sqm7-14)*ZZ+1.0/21*(3*sqm7-7)}
  },
  // This is A'
  {
    {1.0/18*(-sqm7-3)*ZZ*ZZ-1.0/3*ZZ+1.0/3*(sqm7+1),
     1.0/126*(sqm7-21)*ZZ*ZZ+1.0/21*(-sqm7-7)*ZZ+1.0/14*(-sqm7+7),
     1.0/126*(-sqm7+9)*ZZ*ZZ+1.0/14*(-sqm7+1)*ZZ+1.0/42*(sqm7+1)},
    {1.0/63*(5*sqm7-9)*ZZ*ZZ+1.0/21*(6*sqm7-10)*ZZ+1.0/21*(4*sqm7+6),
     1.0/63*(-2*sqm7+21)*ZZ*ZZ+1.0/7*(-sqm7+7)*ZZ+1.0/42*(5*sqm7+7),
     -4.0/63*sqm7*ZZ*ZZ+1.0/42*(-5*sqm7-7)*ZZ+1.0/42*(3*sqm7+7)},
    {1.0/9*(sqm7+3)*ZZ*ZZ+1.0/3*(sqm7+3)*ZZ+2.0/3,
     -4.0/63*sqm7*ZZ*ZZ+1.0/21*(-3*sqm7+7)*ZZ+1.0/21*(4*sqm7+14),
     1.0/126*(11*sqm7-21)*ZZ*ZZ+1.0/21*(3*sqm7-14)*ZZ+1.0/21*sqm7}
  }
};
Mtx_t Gens0[2*MaxNoGens];  // The given generators and their inverses
Mtx_t Gens[2*MaxNoGens];   // The conjugated versions of Gens0

#define EltsDesired 30000 // (Approximately) the number of good points we want
#define EltsSiz (2*(GpKSiz+MaxNoGens)*EltsDesired) // Size of full array of points
Elt_t Elts[EltsSiz]; // Elements as generated
Elt_t SavElts[EltsDesired]; // Elements from previous list

int main(int argc,char *(argv[])) {
  int Ix,Ix2,EltIx,NOldElts,NElts2,NElts,Iter,EqualFg,SavedFg;
  Elt_t NewElt,*Elt_p;
  char letters00[NoGens00+1]="UBA",letters[NoGens],
    Kletters[GpKSiz+1]="1";
                   
  EltsFile=fopen(argv[1],"w");
  assert(EltsFile != NULL);

  hsnorms=fopen("hsnorms","w");
  assert(hsnorms != NULL);

  printf("FF0=\n");  PrintMtx(FF0);
  printf("FF=\n");  PrintMtx(FF);
  printf("\nU0=\n");  PrintMtx(Gens00[0]);
  printf("B0=\n");  PrintMtx(Gens00[1]);
  printf("A0=\n");  PrintMtx(Gens00[2]);

  Ix2=0;
  for(Ix=0;Ix<NoGens00;Ix++) {
    memcpy(Gens0[Ix2],Gens00[Ix],sizeof(Mtx_t));
    letters[Ix2]=letters00[Ix];
    Ix2++;
    InvMtx(Gens0[Ix2],Gens00[Ix]);
    letters[Ix2]=letters00[Ix]+'a'-'A';
    Ix2++;
  }
  assert(Ix2 == NoGens);

  PrintMatrices(Gens0,NoGens,"Gens0");

  for(Ix=0;Ix<NoGens;Ix++)assert(CheckUnitary(Gens0[Ix],FF0));
  
  for(Ix=0;Ix<NoGens;Ix++)
    ConjugateMtx(Gens[Ix],Gens0[Ix]);

  PrintMatrices(Gens,NoGens,"Gens");

  for(Ix=0;Ix<NoGens;Ix++)
    assert(CheckUnitary(Gens[Ix],FF));

  MakeGpK(); // Set up the finite group
  // PrintMatrices(GpK,GpKSiz,"GpK");
  for(Ix=0;Ix<GpKSiz;Ix++)assert(CheckUnitary(GpK[Ix],FF));

  Elts[0].p.s2=0.0;
  Elts[0].p.z[0]=0.0;
  Elts[0].p.z[1]=0.0;
  Elts[0].wLength=0;
  memcpy(&Elts[0].w[0],"",1*sizeof(char));
  Elts[0].w[1]=0;
  memcpy(&Elts[0].m,GpK[0],sizeof(Mtx_t)); // The identity matrix
  NOldElts=1;
  SavedFg=FALSE;

  Iter=0;
  while(TRUE) {
    // Apply non-$K$ generators of $\Gamma$ to the old elements
    NOldElts=MIN(NOldElts,EltsSiz/(NoGens+1));
    NElts=NOldElts;
    for(Ix=0;Ix<NoGens;Ix++) {
      for(EltIx=0;EltIx<NOldElts;EltIx++) {
	ApplyGenToElt(&NewElt,Ix,Elts[EltIx],Gens,letters);
	Elt_p=bsearch(&NewElt,&Elts,
		       (size_t) NOldElts,sizeof(Elt_t),CompareElts);
	if(Elt_p == NULL) {
	  assert(NElts<EltsSiz);
	  Elts[NElts]=NewElt;
	  NElts++;
	}
      }
    }
    // Sort all the points
    qsort(&Elts,(size_t) NElts,sizeof(Elt_t),CompareElts);
    // Eliminate duplicates
    NElts2=1;
    for(EltIx=1;EltIx<NElts;EltIx++)
      if(CompareElts(&Elts[EltIx],&Elts[EltIx-1])>0) {
	Elts[NElts2]=Elts[EltIx];
	NElts2++;
      }
    printf("After Gens: Iter=%i, NOldElts=%i, NElts=%i, NElts2=%i\n",
	   Iter,NOldElts,NElts,NElts2);
    NOldElts=NElts2;

    // Apply elements of $K$ to the old points
    NOldElts=MIN(NOldElts,EltsSiz/(GpKSiz+1));
    NElts=NOldElts;
    for(Ix=0;Ix<GpKSiz;Ix++)
      for(EltIx=0;EltIx<NOldElts;EltIx++) {
	ApplyGenToElt(&NewElt,Ix,Elts[EltIx],GpK,Kletters);
	Elt_p=bsearch(&NewElt,&Elts,
		      (size_t) NOldElts,sizeof(Elt_t),CompareElts);
	if(Elt_p == NULL) {
	  assert(NElts<EltsSiz);
	  Elts[NElts]=NewElt;
	  NElts++;
	}
    }

    // Sort all the elements
    qsort(&Elts,(size_t) NElts,sizeof(Elt_t),CompareEltsWithLength);
    // Eliminate duplicates
    NElts2=1;
    for(EltIx=1;EltIx<NElts;EltIx++)
      if(CompareElts(&Elts[EltIx],&Elts[EltIx-1])>0) {
	Elts[NElts2]=Elts[EltIx];
	NElts2++;
      }
    printf("After  KGp: Iter=%i, NOldElts=%i, NElts=%i, NElts2=%i\n",
	   Iter,NOldElts,NElts,NElts2);
    NOldElts=NElts2;

    // If we have a set of saved elements available, and if our new
    // points match those, we're done.
    if(SavedFg && memcmp(&SavElts,&Elts,sizeof(SavElts)) == 0) break;

    // If we have at least EltsDesired points available, make a copy
    // of them.
    if(NOldElts >= EltsDesired) {
      memcpy(&SavElts,&Elts,sizeof(SavElts));
      SavedFg=TRUE;
    }
    Iter++;
  }
  // printf("\n\n===== * =====\n\n");
  // PrintElts(SavElts,EltsDesired);
  OutputElts(EltsFile,SavElts,EltsDesired);
}