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


void GetGenElts(FILE *file,int *NoElts_p,Elt_t **Elts);

int main(int argc,char *(argv[])) {
  int NoGenElts,GenEltIx,KSiz,EltsSiz,EltsDesired;
  Elt_t *Elts,*GenElts,*KElts,NewElt,*Elt_p,*SavElts;
  int *KMult,*KInv,MaxWordLen;
  int Ix,EltIx,EltIx2,NOldElts,NElts2,NElts,Iter,EqualFg,SavedFg;
  FLOAT_t d2Limit;
  Pt_t Pt;

  // Input files
  FILE *GensFile,*KFile;

  // Output files
  FILE *EltsFile;

  char *Prefix_p,FileNam[MAX_FILENAME_LEN]="../",*Desired_p;

  // Get the prefix
  Prefix_p=getenv("PREFIX");
  assert(Prefix_p != NULL);
  // Get the number of desired elements
  Desired_p=getenv("DESIRED");
  assert(Desired_p != NULL);
  assert(sscanf(Desired_p,"%i",&EltsDesired) == 1);

  // Open the files
  {
    size_t len=strlen(Prefix_p),lenDesired=strlen(Desired_p);
    assert(3+len+1+lenDesired+6<MAX_FILENAME_LEN);
    strcpy(&FileNam[3],Prefix_p);
    // Input files
    strcpy(&FileNam[3]+len,"-Gens");
    printf("%s\n",FileNam);
    GensFile=fopen(FileNam,"r");
    assert(GensFile != NULL);
    strcpy(&FileNam[3]+len,"-K");
    printf("%s\n",FileNam);
    KFile=fopen(FileNam,"r");
    assert(KFile != NULL);
    // Output file
    strcpy(&FileNam[3]+len,"-");
    strcpy(&FileNam[3]+len+1,Desired_p);
    strcpy(&FileNam[3]+len+1+lenDesired,"Elts");
    printf("%s\n",FileNam);
    EltsFile=fopen(FileNam,"w");
    assert(EltsFile != NULL);
  }
  
  // Get the elements of K
  GetKEltsKK(KFile,&KSiz,&KElts,&KMult,&KInv);

  // get the non-K generating elements
  GetGenElts(GensFile,&NoGenElts,&GenElts);
  EltsSiz=4*EltsDesired;
  printf("EltSiz=%i EltsDesired=%i\n",EltsSiz,EltsDesired);
  Elts=malloc(EltsSiz*sizeof(Elt_t));
  assert(Elts != NULL);
  SavElts=malloc(EltsDesired*sizeof(Elt_t));  
  assert(SavElts != NULL);

  // The identity element
  memcpy(&Elts[0],&KElts[0],sizeof(Elt_t));

  d2Limit=HUGE_FLOAT;
  NOldElts=1;
  SavedFg=FALSE;

  // In the first phase we apply generators (in and out of
  // $K$) to the current list of elements.
  Iter=0;
  while(Iter < 50) {
    // Apply non-$K$ generators of $\Gamma$ to the old elements
    NOldElts=MIN(NOldElts,2*EltsDesired);
    NElts=NOldElts;
    for(Ix=0;Ix<NoGenElts;Ix++) {
      for(EltIx=0;EltIx<NOldElts;EltIx++) {
	ApplyEltToEltKK(&NewElt,GenElts[Ix],Elts[EltIx],
			KElts,KMult,KSiz);
	if(NewElt.d2<d2Limit) {
	  Elt_p=bsearch(&NewElt,Elts,
			(size_t) NOldElts,sizeof(Elt_t),CompareElts);
	  if(Elt_p == NULL) {
	    if(NElts >= EltsSiz) break;
	    Elts[NElts]=NewElt;
	    NElts++;
	  }
	  else if(CompareEltsWithWords(&NewElt,Elt_p)<0)
	    memcpy(Elt_p,&NewElt,sizeof(Elt_t));
	}
      }
      if(NElts >= EltsSiz) break;
    }
    // Sort all the points
    qsort(Elts,(size_t) NElts,sizeof(Elt_t),CompareEltsWithWords);
    // 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, d2Limit="FOCnv
	   ", NOldElts=%i, NElts=%i, NElts2=%i\n",
	   Iter,d2Limit,NOldElts,NElts,NElts2);
    fflush(stdout);
    NOldElts=NElts2;

    // Apply elements of $K$ to the old points
    NOldElts=MIN(NOldElts,2*EltsDesired);
    NElts=NOldElts;
    for(EltIx=0;EltIx<NOldElts;EltIx++) {
      for(Ix=1;Ix<KSiz;Ix++) {
	ApplyEltToEltKK(&NewElt,KElts[Ix],Elts[EltIx],
			KElts,KMult,KSiz);
	if(NewElt.d2<d2Limit) {
	  Elt_p=bsearch(&NewElt,Elts,
			(size_t) NOldElts,sizeof(Elt_t),CompareElts);
	  if(Elt_p == NULL) {
	    if(NElts >= EltsSiz) break;	
	    Elts[NElts]=NewElt;
	    NElts++;
	  }
	  else if(CompareEltsWithWords(&NewElt,Elt_p)<0)
	    memcpy(Elt_p,&NewElt,sizeof(Elt_t));
	}
	if(NElts >= EltsSiz) break;
      }
    }
    // Sort all the elements
    qsort(Elts,(size_t) NElts,sizeof(Elt_t),CompareEltsWithWords);
    // 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, d2Limit="FOCnv
	   ", NOldElts=%i, NElts=%i, NElts2=%i\n",
	   Iter,d2Limit,NOldElts,NElts,NElts2);
    fflush(stdout);
    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,EltsDesired*sizeof(Elt_t)) == 0))
      break;

    // If we have at least EltsDesired points available, make a copy
    // of them.
    if(NOldElts >= EltsDesired) {
      memcpy(SavElts,Elts,EltsDesired*sizeof(Elt_t));
      d2Limit=SavElts[EltsDesired-1].d2+epsBigger;
      SavedFg=TRUE;
    }
    Iter++;
  }


  // In the second phase, having obtained (we hope) a reasonably
  // complete list of the first EltsDesired elements, we try to
  // complete the list by multiplying those elements by one another.
  Iter=0;
  while(Iter < 10) {
    NElts=EltsDesired;
    for(EltIx=1;EltIx<EltsDesired;EltIx++) {
      for(EltIx2=0;EltIx2<EltsDesired;EltIx2++) {
	ApplyMtxToPt(Pt,Elts[EltIx].Mtx,Elts[EltIx2].Pt);
	if(PtAbsSq(Pt)<d2Limit) {
	  ApplyEltToEltKK(&NewElt,Elts[EltIx],Elts[EltIx2],
			  KElts,KMult,KSiz);
	  assert(EqPt(NewElt.Pt,Pt,epsSmall));
	  Elt_p=bsearch(&NewElt,Elts,
			(size_t) EltsDesired,sizeof(Elt_t),CompareElts);
	  if(Elt_p == NULL) {
	    if(NElts >= EltsSiz) break;	
	    Elts[NElts]=NewElt;
	    NElts++;
	  }
	  else if(CompareEltsWithWords(&NewElt,Elt_p)<0)
	    memcpy(Elt_p,&NewElt,sizeof(Elt_t));
	}
      }
      if(NElts >= EltsSiz) break;	
    }
    // Sort all the points
    qsort(Elts,(size_t) NElts,sizeof(Elt_t),CompareEltsWithWords);
    // 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  All: Iter=%i, d2Limit="FOCnv
	   ", NElts=%i, NElts2=%i\n",
	   Iter,d2Limit,NElts,NElts2);
    fflush(stdout);
    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,EltsDesired*sizeof(Elt_t)) == 0))
      break;

    // If we have at least EltsDesired points available, make a copy
    // of them.
    memcpy(SavElts,Elts,EltsDesired*sizeof(Elt_t));
    d2Limit=SavElts[EltsDesired-1].d2+epsBigger;
    SavedFg=TRUE;
    Iter++;
  }

  MaxWordLen=0;
  for(EltIx=0;EltIx<EltsDesired;EltIx++)
    if(strlen(SavElts[EltIx].Word)>MaxWordLen)
      MaxWordLen=strlen(SavElts[EltIx].Word);
  printf("Maximum Word Length=%i\n",MaxWordLen);

  fprintf(EltsFile,"NumElts=%i\n",EltsDesired);
  for(EltIx=0;EltIx<EltsDesired;EltIx++) 
    OutputElt(EltsFile,SavElts[EltIx]);
}

void GetGenElts(FILE *ifile,int *NoElts_p,Elt_t **Elts) {
  int EltIx;

  assert(fscanf(ifile,"NoGenElts=%i ",NoElts_p) == 1);
  *Elts=malloc((*NoElts_p)*sizeof(Elt_t));

  for(EltIx=0;EltIx<*NoElts_p;EltIx++) {
    InputElt(ifile,&(*Elts)[EltIx]);
  }
  assert(fscanf(ifile,"%*[ \n\t]") == EOF);
}