#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 finds fixed points of elements of the group

int main(int argc,char *(argv[])) {
  int NumElts,NumFewElts,NumPts,FixerSiz;
  int EltIx1,EltIx2,PtIx,FixerIx,Fixer[MaxFixerSiz];
  int OKFg,EltOrder;

  Pt_t NewPt,*Pts;
  Mtx_t Prod,Inv;
  Elt_t TestElt,*Elt_p;

  // Input files
  FILE *FewEltsFile;

  // Output files

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

  Elt_t *FewElts;
  FiniteOrder_t *FewEltsFO;

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

  // 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");
    // Copy the points onto the points list
  SetupPts(&Pts,&NumPts,FewElts,NumFewElts);
  printf("NumFewElts=%i NumPts=%i\n",NumFewElts,NumPts);
  fflush(stdout);

  // Work out which of the elements are of finite order
  FewEltsFO=malloc(NumFewElts*sizeof(FiniteOrder_t));
  assert(FewEltsFO != NULL);
  for(EltIx1=0;EltIx1<NumFewElts;EltIx1++)
    CalculateFiniteOrder(&FewEltsFO[EltIx1],FewElts[EltIx1]);

  // Find fixed points in the Dirichlet domain
  for(EltIx1=0;EltIx1<NumFewElts;EltIx1++) {
    // Is FewElts[EltIx1] of finite order?
    if(FewEltsFO[EltIx1].EltOrder == 0)continue;
    // Yes. Is the fixed point in the Dirichlet fundamental domain?
    if(! InDirichlet(FewEltsFO[EltIx1].FixedPt,Pts,NumPts))continue;
    // Yes. Find other elements fixing that same fixed point.
    FixerIx=0;
    OKFg=TRUE;
    for(EltIx2=0;EltIx2<NumFewElts;EltIx2++) {
      ApplyMtxToPt(NewPt,FewElts[EltIx2].Mtx,FewEltsFO[EltIx1].FixedPt);
      // Does FewElts[EltIx2] fix this same fixed point?
      if(EqPt(NewPt,FewEltsFO[EltIx1].FixedPt,epsBigger)) {
	// We wish to consider this fixed point now only if (a)
	// FewElts[EltIx1] is the first non-reflection element which
	// fixes it or (b) it is fixed only by complex reflections,
	// and FewElts[EltIx1] is the first of these.

	// Is FewElts[EltIx2] a complex reflection?
	if(FewEltsFO[EltIx2].Phase[0] == 0) {
	  // Yes.  If FewElts[EltIx1] is also a complex relection, and
	  // if FewElts[EltIx2] precedes FewElts[EltIx1] then we don't
	  // proceed.
	  if(FewEltsFO[EltIx1].Phase[0] == 0 && EltIx2<EltIx1)
	    OKFg=FALSE;
	}
	else {
	  // No.  If FewElts[EltIx1] is a complex reflection, we
	  // don't proceed.  If FewElts[EltIx2] precedes
	  // FewElts[EltIx1], we don't precede.
	  
	  //printf("A: EltIx1=%i EltIx2=%i "
	  //	 "FewEltsFO[EltIx1].Phase =%i %i "
	  //	 "FewEltsFO[EltIx2].Phase =%i %i\n",
	  //	 EltIx1,EltIx2,
	  //	 FewEltsFO[EltIx1].Phase[0],FewEltsFO[EltIx1].Phase[1],
	  //	 FewEltsFO[EltIx2].Phase[0],FewEltsFO[EltIx2].Phase[1]);
	  if(FewEltsFO[EltIx1].Phase[0] == 0 || EltIx2 < EltIx1)
	    OKFg=FALSE;
	  //printf("B: OKFg=%i\n",OKFg);
	}
	if(!OKFg)break;
	Fixer[FixerIx]=EltIx2;
	FixerIx++;
      }
      // Does FewElts[EltIx2] shift this fixed point to a point
      // which is still in the fundamental domain, and ``smaller'' than
      // the original (relative to an arbitrary order)?  If so, we'll
      // skip this fixed point, since the same conjugacy class of
      // fixing subgroups is given by the ``smaller'' point.

      if(ComparePts(&NewPt,&(FewEltsFO[EltIx1].FixedPt))<0 && 
	 InDirichlet(NewPt,Pts,NumPts)) {
	
	//InvMtx(Inv,FewElts[EltIx2].Mtx);
	//MultMtx(Prod,FewElts[EltIx2].Mtx,FewElts[EltIx1].Mtx);
	//MultMtx(TestElt.Mtx,Prod,Inv);
	//FixUpElt(&TestElt);
	//Elt_p=bsearch(&TestElt,FewElts,
	//	      (size_t) NumFewElts,sizeof(Elt_t),CompareElts);
	//assert(Elt_p != NULL);

	OKFg=FALSE;
	break;
      }
    }
    if(! OKFg)continue;
    FixerSiz=FixerIx;

    printf("=== Fixed Pt: FixerSiz=%i\n",FixerSiz);
    OutputPt(stdout,FewEltsFO[EltIx1].FixedPt);
    printf("Fixer=\n");
    for(FixerIx=0;FixerIx<FixerSiz;FixerIx++) {
      EltIx2=Fixer[FixerIx];
      printf("%s :   Eigenvalues = e(%i/%i)  e(%i/%i)\n",
	     FewElts[EltIx2].Word,
	     FewEltsFO[EltIx2].Phase[0],FewEltsFO[EltIx2].EltOrder,
	     FewEltsFO[EltIx2].Phase[1],FewEltsFO[EltIx2].EltOrder);
    }
    printf("\n");
  }
}