#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"
#include "../numeric/SU21-fundamental-domain.c"
#include "../numeric/SU21-order.c"

// This program generates a list of conjugates of each element of
// finite order.

int main(int argc,char *(argv[])) {
  int NumFewElts,NumManyElts;
  int FewEltIx,ManyEltIx,EltIx;
  int *ConjFg,OKFg;
  char (*ConjWord)[MaxWordLength+1];

  Elt_t *FewElt_p,*FewElt2_p,*Elt_p,TestElt;
  Mtx_t Prod,Inv;

  // Input files
  FILE *FewEltsFile;
  FILE *EltsFile;
  FILE *ManyEltsFile;

  // Output files

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

  Elt_t *FewElts,*ManyElts;
  FiniteOrder_t FO;

  // Get the prefix
  Prefix_p=getenv("PREFIX");
  assert(Prefix_p != NULL);

  // Open the files
  {
    size_t len=strlen(Prefix_p);
    assert(3+len+9<MAX_FILENAME_LEN);
    strcpy(&FileNam[3],Prefix_p);
    // Input files

    strcpy(&FileNam[3]+len,"-FewElts");
    printf("%s\n",FileNam);
    FewEltsFile=fopen(FileNam,"r");
    assert(FewEltsFile != NULL);
    fflush(stdout);

    strcpy(&FileNam[3]+len,"-ManyElts");
    printf("%s\n",FileNam);
    ManyEltsFile=fopen(FileNam,"r");
    assert(ManyEltsFile != NULL);
    fflush(stdout);
  }

  // Get the few elements
  GetEltsString(FewEltsFile,"NumFewElts",&NumFewElts,&FewElts);
  qsort(FewElts,(size_t) NumFewElts,sizeof(Elt_t),CompareEltsWithWords);
  if(strcmp(FewElts[0].Word,"") == 0)
    strcpy(FewElts[0].Word,"1");
  // Get the many elements
  GetEltsString(ManyEltsFile,"NumElts",&NumManyElts,&ManyElts);
  qsort(ManyElts,(size_t) NumManyElts,sizeof(Elt_t),CompareEltsWithWords);
  if(strcmp(ManyElts[0].Word,"") == 0)
    strcpy(ManyElts[0].Word,"1");

  ConjFg=malloc(NumManyElts*sizeof(int));
  assert(ConjFg != NULL);
  ConjWord=malloc(NumManyElts*(MaxWordLength+1)*sizeof(char));
  assert(ConjWord != NULL);

  FewElt_p=FewElts;
  for(FewEltIx=0;FewEltIx<NumFewElts;FewEltIx++,FewElt_p++) {
    // Deal with *FewElt_p only if it has finite order
    if(Order(FewElt_p->Mtx) == 0) continue;
    CalculateFiniteOrder(&FO,*FewElt_p);
    // The ConjFg array will be set TRUE at positions corresponding to
    // conjugates of *FewElt_p
    memset(ConjFg,FALSE,NumManyElts*sizeof(int));
    OKFg=TRUE;
    for(ManyEltIx=0;ManyEltIx<NumManyElts;ManyEltIx++) {
      // See what happens when we conjugate by ManyElts[ManyEltIx]
      InvMtx(Inv,ManyElts[ManyEltIx].Mtx);
      MultMtx(Prod,Inv,FewElt_p->Mtx);
      MultMtx(TestElt.Mtx,Prod,ManyElts[ManyEltIx].Mtx);
      FixUpElt(&TestElt);
      // If the new conjugate has already been studied, there's no
      // point in redoing this conjugacy class.
      FewElt2_p=bsearch(&TestElt,FewElts,
			(size_t) NumFewElts,sizeof(Elt_t),CompareElts);
      if(FewElt2_p != NULL && (FewElt_p-FewElt2_p) > 0) {
	OKFg=FALSE;
	break;
      }
      // Look for the new conjugate in ManyElts
      Elt_p=bsearch(&TestElt,ManyElts,
		    (size_t) NumManyElts,sizeof(Elt_t),CompareElts);
      // If it's there, mark it as a conjugate and save the
      // conjugating element.
      if(Elt_p != NULL && !ConjFg[Elt_p-ManyElts]) {
	ConjFg[Elt_p-ManyElts]=TRUE;
	strcpy(ConjWord[Elt_p-ManyElts],ManyElts[ManyEltIx].Word);
      }
    }
    if(!OKFg)continue;

    printf("=== %s ;  Eigenvalues = e(%i/%i)  e(%i/%i)\n",
	   FewElt_p->Word,
	   FO.Phase[0],FO.EltOrder,
	   FO.Phase[1],FO.EltOrder);
    for(EltIx=0;EltIx<NumManyElts;EltIx++)
      if(ConjFg[EltIx])
	printf("%s:(%s)^(%s)\n",
	       ManyElts[EltIx].Word,FewElt_p->Word,ConjWord[EltIx]);
    printf("\n");
  }
}