// Calling sequence
//   
//  PREFIX=cc/cc-dd ./permsearch
//

#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 "../permsearch/SU21-permsearch-base.c"

// This program has to do with permutations~$B$ which act on the left
// on a group~$K$ and satisfy various relations of the form:
//
//    B k_1 = k_2 B  and B k_3 = k_4 B^-1
//
// for various given pairs of elements~$(k_1,k_2)$ and $(k_3,k_4)$
// of~$K$, where the~$k_i$ act on~$K$ by left-multiplication.
//
// The set on which $B$~acts, which is a copy of the group~$K$, enters
// into the conditions only through its structure as a left-$K$-set.
// Consequently any permutation of~$K$ which preserves its structure
// as a left $K$-set will conjugate any valid~$B$ to another such.  It
// is obvious (and very well-known) that a permutation preserves the
// structure of~$K$ as a left-$K$-set if and only if it is
// right-multiplication by some~$k\in K$.
//
// If such a B satisfies B(1)=y, then it will also satisfy
//
//   B(k1)    = B k1 (1) = k2 B (1)    = k2 y
//
//   B(k3 y)  = B k3 (y) = k4 B^-1 (y) = k4 1 = k4
//
// If B is conjugated by right-translation by k1 (respectively k3 y)
// this gives equivalent actions satisfying:
//
//   B'(1)  = k2 y k1^-1
//
//   B''(1) = k4 y^-1 k3^-1
//
// So, if one is only interested in one action from each equivalence class,
// one may consider only B's for which B(1)=y, and skip B's for which
//
//   B(1)= k2 y k1^-1    or    B(1) = k4 y^-1 k3^-1
//
// This program prints a list of possibilities for such shortcuts.

int CompareInts(const void *Ptr1,const void *Ptr2);

int main(int argc,char *(argv[])) {
  // Information about the group $K$
  int KSiz,*KMult,*KInv;
  Elt_t *KElts;
  // Information about the relations
  int NumRels;
  Rels_t Rels;

  // Input files
  FILE *KFile;
  FILE *KRelFile;
  FILE *KGAPFile;

  INT_t *k1Perm,*k2Perm,*k3Perm,*k4Perm;

  // List of ``equivalent'' y's
  int *yEquiv,x,y,RelIx,Ix;

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

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

  // Open the files
  {
    size_t len=strlen(Prefix_p);
    assert(3+len+7<MAX_FILENAME_LEN);
    strcpy(&FileNam[3],Prefix_p);
    // Input files
    strcpy(&FileNam[3]+len,"-K");
    printf("%s\n",FileNam);
    KFile=fopen(FileNam,"r");
    assert(KFile != NULL);
    strcpy(&FileNam[3]+len,"-KRel");
    printf("%s\n",FileNam);
    KRelFile=fopen(FileNam,"r");
    assert(KRelFile != NULL);
    // Output files
    strcpy(&FileNam[3]+len,"-K-gap");
    printf("%s\n",FileNam);
    KGAPFile=fopen(FileNam,"w");
    assert(KGAPFile != NULL);
    fflush(stdout);
  }

  // Get the elements of~$K$
  GetKEltsKK(KFile,&KSiz,&KElts,&KMult,&KInv);
  assert(strcmp(KElts[0].Word,"") == 0 ||
	 strcmp(KElts[0].Word,"1") == 0);
  strcpy(KElts[0].Word,"1");
  assert(KSiz < SHRT_MAX);
  // Set up a file to check the K just found using GAP
  fprintf(KGAPFile,"KList:=[\n");
  x=0;
  fprintf(KGAPFile,"(%-s)",KElts[x].Word);
  for(x=1;x<KSiz;x++) 
    fprintf(KGAPFile,",\n(%-s)",KElts[x].Word);
  fprintf(KGAPFile,"\n];;\n\n");

  fprintf(KGAPFile,"KMultList:=[\n");
  for(x=0;x<KSiz;x++) {
    for(y=0;y<KSiz;y++) {
      fprintf(KGAPFile,"[(%s),(%s),(%s)]",
	      KElts[x].Word,KElts[y].Word,
	      KElts[KMult[KSiz*x+y]].Word);
      if(y != KSiz-1 || x != KSiz-1)
	fprintf(KGAPFile,",\n");
      else
	fprintf(KGAPFile,"\n];;\n\n");
    }
    if(x != KSiz-1)fprintf(KGAPFile,"\n");
  }
  
  fprintf(KGAPFile,"KInvList:=[\n");
  x=0;
  fprintf(KGAPFile,"[(%s),(%s)]",
	  KElts[x].Word,KElts[KInv[x]].Word);
  for(x=1;x<KSiz;x++)
    fprintf(KGAPFile,",\n[(%s),(%s)]",
	    KElts[x].Word,KElts[KInv[x]].Word);
  fprintf(KGAPFile,"\n];;\n");


  // Get the pairs~$(k,k')$ which constrain~$B$
  GetKRels(KRelFile,&NumRels,&Rels,KSiz,KElts,KMult);

  yEquiv=malloc(2*NumRels*sizeof(int));
  assert(yEquiv != NULL);

  for(y=0;y<KSiz;y++) {
    // Find a set of equivalent y's
    k1Perm=Rels.k1; k2Perm=Rels.k2; k3Perm=Rels.k3; k4Perm=Rels.k4;
    Ix=0;
    for(RelIx=0;RelIx<NumRels;RelIx++) {
      // B(1)= k2 y k1^-1
      yEquiv[Ix]=
	KMult[KSiz*k2Perm[0]+KMult[KSiz*y+KInv[k1Perm[0]]]];
      Ix++;
      // B(1) = k4 y^-1 k3^-1
      yEquiv[Ix]=
	KMult[KSiz*k4Perm[0]+KMult[KSiz*KInv[y]+KInv[k3Perm[0]]]];
      Ix++;

      k1Perm+=KSiz; k2Perm+=KSiz; k3Perm+=KSiz; k4Perm+=KSiz; 
    }
    // Sort them
    qsort(yEquiv,(size_t) 2*NumRels,sizeof(int),CompareInts);
    // Print them
    for(Ix=0;Ix<2*NumRels;Ix++)
      printf("%4i",yEquiv[Ix]);
    printf("\n");
  }
}

// Compare two integers
int CompareInts(const void *Ptr1,const void *Ptr2) {
  const int *Int1=(const int *)Ptr1, *Int2=(const int *)Ptr2;

  if(*Int1<*Int2)return -1;
  if(*Int1>*Int2)return +1;
  return 0;
}