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


// Here we try to find representatives of the first so many
// K-double-cosets of elements of our group.  Here K is supposed to be
// the subgroup fixing the origin.

void GetGenElts(FILE *file,int *NoElts_p,Elt_t **Elts);
inline void KOrbitRep(Pt_t Rep,int *KIx,Pt_t Pt,
		      Elt_t KElts[],int KSiz);

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

  // Input files
  FILE *GensFile,*KFile;

  // Output files
  FILE *CosetsFile;

  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+11<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,"KKCosets");
    printf("%s\n",FileNam);
    CosetsFile=fopen(FileNam,"w");
    assert(CosetsFile != NULL);
  }
  
  // Get the elements of K
  GetKEltsKK(KFile,&KSiz,&KElts,&KMult,&KInv);
  // Check that they fix the origin
  for(Ix=0;Ix<KSiz;Ix++) {
    ApplyMtxToPt(Pt,KElts[Ix].Mtx,ZeroPt);
    assert(sqrt(PtAbsSq(Pt))<epsSmall);
  }

  // get the non-K generating elements
  GetGenElts(GensFile,&NoGenElts,&GenElts);
  
  // Allocate memory for the elements list
  EltsSiz=(NoGenElts+1)*2*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
  NOldElts=1;
  memcpy(&Elts[0],&KElts[0],sizeof(Elt_t));

  // In the first phase we repeatedly apply generators to our K-double
  // cosets and try to get a reasonable approximation to the desired
  // list.
  d2Limit=HUGE_FLOAT;
  SavedFg=FALSE;

  Iter=0;
  while(Iter < 50) {
    // Apply non-$K$ generators of $\Gamma$ to the old elements
    NOldElts=MIN(NOldElts,2*EltsDesired);
    // if(NOldElts == 2*EltsDesired)
    //   d2Limit=Elts[2*EltsDesired-1].d2+epsBigger;
    NElts=NOldElts;
    for(Ix=0;Ix<NoGenElts;Ix++) {
      for(EltIx=0;EltIx<NOldElts;EltIx++) {
	ApplyMtxToPt(Pt,GenElts[Ix].Mtx,Elts[EltIx].Pt);
	d2=PtAbsSq(Pt);
	if(d2<d2Limit) {
	  KOrbitRep(Rep,&KIx,Pt,KElts,KSiz);
	  ApplyEltToEltKK(&tmpElt,GenElts[Ix],Elts[EltIx],
			  KElts,KMult,KSiz);
	  ApplyEltToEltKK(&NewElt,KElts[KIx],tmpElt,
			  KElts,KMult,KSiz);
	  assert(EqPt(NewElt.Pt,Rep,epsSmall));
	  Elt_p=bsearch(&NewElt,Elts,(size_t) NOldElts,sizeof(Elt_t),
			CompareEltsPtsOnly);
	  if(Elt_p == NULL) {
	    assert(NElts < EltsSiz);
	    Elts[NElts]=NewElt;
	    NElts++;
	  }
	  else if(CompareEltsWithWords(&NewElt,Elt_p)<0)
	    memcpy(Elt_p,&NewElt,sizeof(Elt_t));
	}
      }
    }
    // Sort all the elements
    qsort(Elts,(size_t) NElts,sizeof(Elt_t),CompareEltsWithWords);
    // Eliminate duplicates
    NElts2=1;
    for(EltIx=1;EltIx<NElts;EltIx++)
      if(CompareEltsPtsOnly(&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;

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

  //fprintf(CosetsFile,"FIRST PHASE ONLY: NumManyKKCosets=%i\n",EltsDesired);
  //for(EltIx=0;EltIx<EltsDesired;EltIx++) 
  //  OutputElt(CosetsFile,SavElts[EltIx]);


  // In the second phase, having obtained (we hope) a reasonably
  // complete list of K-double-coset representatives, we complete the
  // list by multiplying them by one another
  d2Limit=Elts[EltsDesired-1].d2+epsBigger;
  Iter=0;
  while(Iter < 10) {
    NElts=EltsDesired;
    for(EltIx=1;EltIx<EltsDesired;EltIx++) {
      for(KIx1=0;KIx1<KSiz;KIx1++) {
	ApplyMtxToPt(Pt1,KElts[KIx1].Mtx,Elts[EltIx].Pt);
	for(EltIx2=0;EltIx2<EltsDesired;EltIx2++) {
	  ApplyMtxToPt(Pt,Elts[EltIx2].Mtx,Pt1);
	  d2=PtAbsSq(Pt);
	  if(d2<d2Limit) {
	    KOrbitRep(Rep,&KIx,Pt,KElts,KSiz);
	    ApplyEltToEltKK(&NewElt,KElts[KIx1],Elts[EltIx],
			    KElts,KMult,KSiz);
	    ApplyEltToEltKK(&tmpElt,Elts[EltIx2],NewElt,
			    KElts,KMult,KSiz);
	    ApplyEltToEltKK(&NewElt,KElts[KIx],tmpElt,
			    KElts,KMult,KSiz);
	    //assert(EqPt(NewElt.Pt,Rep,epsBiggest));

	    Elt_p=bsearch(&NewElt,Elts,(size_t) EltsDesired,sizeof(Elt_t),
			  CompareEltsPtsOnly);
	    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;
      }
      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(CompareEltsPtsOnly(&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(CosetsFile,"NumKKCosets=%i\n",EltsDesired);
  for(EltIx=0;EltIx<EltsDesired;EltIx++) 
    OutputElt(CosetsFile,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);
}

// Find the ``smallest'' point in the K-orbit of a given point.
inline void KOrbitRep(Pt_t Rep,int *KIx,Pt_t Pt,
		      Elt_t KElts[],int KSiz) {
  int Ix;
  Pt_t TrialPt;
  memcpy(Rep,Pt,sizeof(Pt_t));
  *KIx=0;
  for(Ix=1;Ix<KSiz;Ix++) {
    ApplyMtxToPt(TrialPt,KElts[Ix].Mtx,Pt);
    if(ComparePts(TrialPt,Rep) < 0) {
      //printf("C:\n");
      *KIx=Ix;
      memcpy(Rep,TrialPt,sizeof(Pt_t));
    }
  }
  //if(ComparePts(TrialPt,Pt)<0) {
  //  printf("Pt =");
  //  OutputPt(stdout,Pt);
  //  printf("Rep=");
  //  OutputPt(stdout,Rep);
  //}
}