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

#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,EltIx,NOldElts,NElts2,NElts,Iter,EqualFg,SavedFg;
  Elt_t NewElt,*Elt_p;
  char letters[4]="UBA",Kletters[2]="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(Gens0[0]);
  printf("A0=\n");  PrintMtx(Gens0[1]);
  printf("B0=\n");  PrintMtx(Gens0[2]);

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