#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <complex.h>
#include <string.h>

#include "a7p2-0-base.c"

typedef struct {
  Pt_t p; // Fixed point
  Elt_t e; // Element which fixes it
  int order;
} FixPair_t;

#define MaxNumElts 30000
Elt_t Elts[MaxNumElts];

void COutputMtx(FILE *MtxFile,Mtx_t mtx);
int CompareFixPairsPtsOnly(const void *Ptr1,const void *Ptr2);
int CompareFixPairs(const void *Ptr1,const void *Ptr2);

int main(int argc,char *(argv[])) {
  int NumElts,Ix,DirIx,DirIx2,NDir,order,NDir2;
  int FixIx,FixIx2,NFix,NFix2,NumPts;
  Pt_t FixedPt,NewPt;
  FILE *MtxFile;
#define MaxNoFixed 1000 // All fixed points
  FixPair_t FixPairs[MaxNoFixed];
#define MaxNoDir 100 // Fixed points in the Dirichlet domain
  FixPair_t DirPairs[MaxNoDir];

  EltsFile=fopen(argv[1],"r");
  assert(EltsFile != NULL);

  MtxFile=fopen("tmpMtx","w");
  assert(MtxFile != NULL);

  NumElts=ReadElts(EltsFile,Elts,MaxNumElts);

  NumPts=MIN(MaxNumPts,NumElts);
  for(Ix=0;Ix<NumPts;Ix++)
    memcpy(&Pts[Ix],&Elts[Ix].p,sizeof(Pt_t));

  // NOTA BENE: The whole approach may fail for complex reflections.
  // It is set up to work when the finite order elements have unique
  // fixed points.  
 
  // Find fixed points
  NFix=NDir=0;
  for(Ix=0;Ix<NumElts;Ix++) {
    order=Order(Elts[Ix].m,&FixedPt);
    if(order > 0) {
      assert(NFix<MaxNoFixed);
      memcpy(&(FixPairs[NFix].p),&FixedPt,sizeof(Pt_t));
      memcpy(&(FixPairs[NFix].e),&Elts[Ix],sizeof(Elt_t));
      FixPairs[NFix].order=order;
      NFix++;
      // See if fixed point is in Dirichlet domain
      if(InDirichlet(FixedPt,NumPts)) {
	assert(NDir<MaxNoDir);
	memcpy(&(DirPairs[NDir].p),&FixedPt,sizeof(Pt_t));
	memcpy(&(DirPairs[NDir].e),&Elts[Ix],sizeof(Elt_t));
        DirPairs[NDir].order=order;
	NDir++;
      }
    }
   }

  // Sort the fixed-point pairs
  qsort(&FixPairs,(size_t) NFix,sizeof(FixPair_t),CompareFixPairs);
  qsort(&DirPairs,(size_t) NDir,sizeof(FixPair_t),CompareFixPairs);


  printf("***** DIRICHLET DOMAIN FINITE ORDER ELEMENTS *****\n");
  for(DirIx=0;DirIx<NDir;DirIx++) {
    printf("\n");
    printf("order = %i\n",DirPairs[DirIx].order);
    OutputPt(stdout,DirPairs[DirIx].p);
    OutputElt(stdout,DirPairs[DirIx].e);
  }

  // Eliminate elements, keeping just one for each fixed points.

  NDir2=1;
  for(DirIx=1;DirIx<NDir;DirIx++)
    if(ComparePts(&DirPairs[DirIx].p,&DirPairs[DirIx-1].p)>0) {
      DirPairs[NDir2]=DirPairs[DirIx];
      NDir2++;
    }
  NDir=NDir2;

  printf("\n***** DIRICHLET DOMAIN FIXED POINTS *****\n");
  for(DirIx=0;DirIx<NDir2;DirIx++) {
    printf("\n");
    printf("order = %i\n",DirPairs[DirIx].order);
    OutputPt(stdout,DirPairs[DirIx].p);
    OutputElt(stdout,DirPairs[DirIx].e);
  }

  // Check whether any of these fixed points lie in common orbits; in
  // other words, check whether the associated finite subgroups are
  // conjugate.
  printf("\n***** DIRICHLET DOMAIN POINTS WHICH ARE CONJUGATE *****\n");
  for(DirIx=0;DirIx<NDir;DirIx++)
    for(DirIx2=0;DirIx2<NDir;DirIx2++)
      if(DirIx2 != DirIx) {
	for(Ix=0;Ix<NumElts;Ix++) {
	  ApplyMtxToPt(&NewPt,Elts[Ix].m,DirPairs[DirIx].p);
	  if(EqPt(NewPt,DirPairs[DirIx2].p)) {
	    printf("\n");
	    OutputPt(stdout,DirPairs[DirIx].p);
	    OutputElt(stdout,DirPairs[DirIx].e);
	    printf(" is in the same orbit as\n");
	    OutputPt(stdout,DirPairs[DirIx2].p);
	    OutputElt(stdout,DirPairs[DirIx2].e);
	    printf("conjugated by\n");
	    OutputElt(stdout,Elts[Ix]);
	    break;
	  }
	}
      }

  printf("\n***** FINITE ORDER ELEMENTS *****\n");
  for(FixIx=0;FixIx<NFix;FixIx++) {
    printf("\n");
    printf("order = %i\n",FixPairs[FixIx].order);
    OutputPt(stdout,FixPairs[FixIx].p);
    OutputElt(stdout,FixPairs[FixIx].e);
  }

  // Eliminate elements, keeping just one for each fixed points.

  NFix2=1;
  for(FixIx=1;FixIx<NFix;FixIx++)
    if(ComparePts(&FixPairs[FixIx].p,&FixPairs[FixIx-1].p)>0) {
      FixPairs[NFix2]=FixPairs[FixIx];
      NFix2++;
    }
  NFix=NFix2;

  printf("\n***** FIXED POINTS *****\n");
  for(FixIx=0;FixIx<NFix2;FixIx++) {
    printf("\n");
    printf("order = %i\n",FixPairs[FixIx].order);
    OutputPt(stdout,FixPairs[FixIx].p);
    OutputElt(stdout,FixPairs[FixIx].e);
  }

  // Check whether any of these fixed points lie in common orbits; in
  // other words, check whether the associated finite subgroups are
  // conjugate.
  printf("\n***** POINTS WHICH ARE CONJUGATE *****\n");
  for(FixIx=0;FixIx<NFix;FixIx++)
    for(FixIx2=0;FixIx2<NFix;FixIx2++)
      if(FixIx2 != FixIx) {
	for(Ix=0;Ix<NumElts;Ix++) {
	  ApplyMtxToPt(&NewPt,Elts[Ix].m,FixPairs[FixIx].p);
	  if(EqPt(NewPt,FixPairs[FixIx2].p)) {
	    printf("\n");
	    OutputPt(stdout,FixPairs[FixIx].p);
	    OutputElt(stdout,FixPairs[FixIx].e);
	    printf(" is in the same orbit as\n");
	    OutputPt(stdout,FixPairs[FixIx2].p);
	    OutputElt(stdout,FixPairs[FixIx2].e);
	    printf("conjugated by\n");
	    OutputElt(stdout,Elts[Ix]);
	    break;
	  }
	}
      }
}

void COutputMtx(FILE *MtxFile,Mtx_t mtx) {
  int ir,ic;

  for(ir=0;ir<3;ir++) {
    if(ir==0)
      fprintf(MtxFile,"  {{");
    else
      fprintf(MtxFile,"   {");
    for(ic=0;ic<3;ic++) {
      if(ic != 0)fprintf(MtxFile,"    ");
      fprintf(MtxFile,"% .15f%+.15f*I",
	      creal(mtx[ir][ic]),cimag(mtx[ir][ic]));
      if(ic != 2)fprintf(MtxFile,",\n");
    }
    if(ir != 2)
      fprintf(MtxFile,"},\n\n");
    else
      fprintf(MtxFile,"}\n  },\n");
  }
  return;
}
	
// Compare two fixed-point pairs, considering, first the fixed points
// themselves, second the orders of the fixing elements, third the
// lengths of the words representing the fixing elements, and finally
// the words themselves
int CompareFixPairs(const void *Ptr1,const void *Ptr2) {
  const FixPair_t *Pair1=(const FixPair_t *)Ptr1,
    *Pair2=(const FixPair_t *)Ptr2;
  int tmp;
  
  tmp=ComparePts(&(Pair1->p),&(Pair2->p));
  if(tmp != 0)return tmp;

  if(Pair1->order > Pair2->order)return +1;
  if(Pair2->order > Pair1->order)return -1;

  if(Pair1->e.wLength > Pair2->e.wLength)return +1;
  if(Pair2->e.wLength > Pair1->e.wLength)return -1;

  return strcmp(Pair1->e.w,Pair2->e.w);
}