/* st06.c CJP Version 1.1 08/10/07 */
#include<stdlib.h>
#include<stdio.h>
#include"blimp.h"

#define MAXGENS 10000

int test=0;
int c;
int *lrel;
int **exrl,**inrl,*extp,*intp,*inv;
int **rellist,*rels,**temp,*tempstore;
int exrels,inrels,ingens;

int comparerels(const void* a, const void* b);
int perm(int* from, int* to);

int main(int argc, char **argv)
{
    int from[MAXGENS];
    int to[MAXGENS];
    int i,j,k;
    int good;
    FILE *f1,*f2;
    int *gp;
    int leave;
    int *pt;
    
    good=1;
    gp = (int*) malloc(sizeof(int)*GRPSIZ);
    if(gp==NULL) good=0;
    intp = (int*) malloc(sizeof(int)*GRPSIZ);
    if(intp==NULL) good=0;
    inv = (int*) malloc(sizeof(int)*GRPSIZ);
    if(inv==NULL) good=0;
    extp = (int*) malloc(sizeof(int)*GRPSIZ);
    if(extp==NULL) good=0;
    exrl = (int**) malloc(sizeof(int*)*GRPSIZ);
    if(exrl==NULL) good=0;
    inrl = (int**) malloc(sizeof(int*)*GRPSIZ);
    if(inrl==NULL) good=0;
    rels = (int*) malloc(sizeof(int)*GRPSIZ);
    if(rels==NULL) good=0;
    rellist = (int**) malloc(sizeof(int*)*GRPSIZ);
    if(rellist==NULL) good=0;
    temp = (int**) malloc(sizeof(int*)*20);
    if(temp==NULL) good=0;
    tempstore = (int*) malloc(sizeof(int)*20);
    if(tempstore==NULL) good=0;
    lrel = (int*) malloc(sizeof(int)*5);
    if(lrel==NULL) good=0;


    if(good==0)
    {
        printf("Error allocating memory in st06\n");
        return 0;
    }
    
    if(argc!=4)
    {
        printf("Usage: st06 <-c or -i> <input> <output>\n");
        return 0;
    }

    f1 = (FILE*) fopen(argv[2],"r");
    if(f1==NULL)
    {
        printf("Error opening file %s\n",argv[2]);
        return 0;
    }
    
    f2 = (FILE*) fopen(argv[3],"w");
    if(f2==NULL)
    {
        printf("Error opening file %s\n",argv[3]);
        return 0;
    }

    if(argv[1][0]=='-' && argv[1][1]=='c')
    {
        c=1;
    }
    else if(argv[1][0]=='-' && argv[1][1]=='i')
    {
        c=0;
    }
    else
    {
        printf("Invalid switch given to st06, use -c or -i\n");
        return 0;
    }
    
    while(1)
    {        
        i = ReadGroup(gp,f1);
        if(i==-1) break;
        
        leave=0;
        // Get Internal and External Presentations
        ExtractBlimp(gp,2,extp);
        ExtractBlimp(gp,6,intp);
        // Make a list of the triangles in the external pres
        exrels=0;
        pt = &extp[4];
         
        for(i=0;i<extp[3];i++)
        {
            j=0;
            while(1)
            {
                if((*pt)!=0)
                {
                    j++;
                    pt++;
                }
                else
                {
                    pt++;
                    break;
                }
                if(j==3 && (*pt)==0)
                {
                    exrl[exrels] = pt-3;
                    exrels++;
                }
            }         
        }
        
        // Sort them
        qsort(exrl,exrels,sizeof(int*),comparerels);
        // Make a list of triangles in internal pres
        inrels=0;
        pt = &intp[3];
        
        for(i=0;i<intp[2];i++)
        {
            j=0;
            while(1)
            {
                if((*pt)!=0)
                {
                    j++;
                    pt++;
                }
                else
                {
                    pt++;
                    break;
                }
                if(j==3 && (*pt)==0)
                {
                    inrl[inrels] = pt-3;
                    inrels++;
                }
            }         
        }
        // If internal has fewer gens, throw away group
        ExtractBlimp(gp,3,inv);
        if(inv[2]<extp[2])
        {
            printf("st06 found fewer internal generators than external ones...\n");
            exit(45);
        }
        // Construct Inverse Table
        k=inv[1];
        ingens=inv[2];
        for(i=0;i<k;i++)
        {
            inv[i] = inv[i+2]; 
        }
        
        // Set Largest Relator
        if(c==1)
        {
            lrel[0]=extp[2];
            lrel[1]=extp[2];
            lrel[2]=extp[2];
        }
        else
        {
            if(exrels>0)
            {
                lrel[0]=exrl[exrels-1][0];
                lrel[1]=exrl[exrels-1][1];
                lrel[2]=exrl[exrels-1][2];
            }
            else
            {
                lrel[0]=extp[2];
                lrel[1]=extp[2];
                lrel[2]=extp[2];
            }
        }
       
        // Do all Permutations, if fail, leave=1 and exit quickly
        for(i=1;i<=ingens;i++)
        {
            from[i]=0;
            to[i]=0;
        }
        
        if(perm(from, to)==1) leave=1;
        // Write out group if no earlier permutation found
        if(leave==1) continue;
        
        WriteGroup(gp,f2);
        
    }
    WriteEOF(f2);
    fclose(f1);
    fclose(f2);    
}

int perm(int* from2, int* to2) // Return 1 reject gp, return 0 reject perm
{
    
    int invol;
    int num,prenum;
    int i,j,k;
    int from[MAXGENS], to[MAXGENS];
    int *pt;
    
    for(i=1;i<=ingens;i++)
    {
        from[i] = from2[i];
        to[i] = to2[i];
    }
    
    // Make list of images of relators
    num=0;
    pt = rels;
    for(i=0;i<inrels;i++)
    {
        if(to[inrl[i][0]]!= 0 && to[inrl[i][1]]!= 0 && to[inrl[i][2]]!= 0)
        {
            rellist[num]=pt;
            *(pt++)=to[inrl[i][0]];
            *(pt++)=to[inrl[i][1]];
            *(pt++)=to[inrl[i][2]];
            num++;
        }
    }
    prenum=num;
    for(i=0;i<6;i++)
    {
        temp[i] = (tempstore + (3*i));
    }
   
    // For each relator in this list
    for(j=0;j<num;j++)
    {
        // Create all rotations
        for(i=0;i<3;i++)
        {
            temp[0][i]=rellist[j][i];
            temp[1][i]=rellist[j][(i+1)%3];
            temp[2][i]=rellist[j][(i+2)%3];
        }
        // Compare them (Finding smallest)
        k=0;
        for(i=0;i<3;i++)
        {
            if(comparerels(&(temp[k]),&(temp[i]))>0) k=i;
        }

        // Write out the smallest
        rellist[j][0]=temp[k][0];
        rellist[j][1]=temp[k][1];
        rellist[j][2]=temp[k][2];
    }
    
    // Sort List  
    qsort(rellist,num,sizeof(int*),comparerels);
    // Remove any higher than largest
    for(i=0;i<num;i++) 
    { 
        if(comparerels(&lrel,&(rellist[i]))<0) 
        { 
            num=i;  
            break; 
        }
    }
    
    // Compare with external list    
    for(i=0;i<num && i<exrels;i++)
    {
        k=comparerels(&rellist[i],&exrl[i]);     
        if(k==0) continue; 
        if(k<0) return 1; 
        if(k>0) return 0;
    }
    
    if(i<num) return 1;
    
    if(i<exrels) // Need to know that we MIGHT make the extra one, else reject perm
    { 
        if(from[exrl[i][0]]!=0 && from[exrl[i][1]]!=0 && from[exrl[i][2]]!=0)
        { 
            return 0;
        } 
    }
    
    if(i==num && i==exrels) // Stop if we're not going to make more rels
    { 
        if(prenum == exrels) return 0;
    } 
    
    // If accepted (fn not returned), Increment Permutation
    // ONLY DO PERMUTATIONS THAT TAKE INVOLUTIONS TO INVOLUTIONS!
    invol=0;
    i=1; 
    while(from[i]) i++;
    if(i==inv[i]) invol=1;
    if(invol==0)
    {
        for(j=1;j<=ingens;j++)
        {
            if(to[j]==0 && j!=inv[j])
            {
                to[j]=i;
                from[i]=j;
                to[inv[j]]=inv[i];
                from[inv[i]]=inv[j];
                
                if(perm(from,to)==1) return 1;
                
                to[j]=0;
                from[i]=0;
                to[inv[j]]=0;
                from[inv[i]]=0;
            }
        }
    }
    if(invol==1)
    {
        for(j=1;j<=ingens;j++)
        {
            if(to[j]==0 && j==inv[j])
            {
                to[j]=i;
                from[i]=j;
                to[inv[j]]=inv[i];
                from[inv[i]]=inv[j];
                
                if(perm(from,to)==1) return 1;
                
                to[j]=0;
                from[i]=0;
                to[inv[j]]=0;
                from[inv[i]]=0;
            }
        }  

    }
    
}


int comparerels(const void *a, const void *b) // Return -ve if first is less
{
    int* zero = *((int**) a);
    int* one = *((int**) b);
    int z0[3],z1[3];
    int wz0,wz1;
    int i,j,k;
    
    for(i=0;i<3;i++)
    {
        if(zero[i]<=zero[(i+1)%3] && zero[i]<=zero[(i+2)%3]) {
            z0[0] = zero[i]; wz0=i; }
        if(one[i]<=one[(i+1)%3] && one[i]<=one[(i+2)%3]) {
            z1[0] = one[i]; wz1 = i; }
    }

    if(zero[(wz0+1)%3] <= zero[(wz0+2)%3])
    {
        z0[1] = zero[(wz0+1)%3];
        z0[2] = zero[(wz0+2)%3];
    }
    else
    {
        z0[1] = zero[(wz0+2)%3];
        z0[2] = zero[(wz0+1)%3];
    }
    if(one[(wz1+1)%3] <= one[(wz1+2)%3])
    {
        z1[1] = one[(wz1+1)%3];
        z1[2] = one[(wz1+2)%3];
    }
    else
    {
        z1[1] = one[(wz1+2)%3];
        z1[2] = one[(wz1+1)%3];
    }
    
    if(z0[2] < z1[2]) return -1;
    if(z0[2] > z1[2]) return 1;
    
    if(z0[1] < z1[1]) return -1;
    if(z0[1] > z1[1]) return 1;

    if(z0[0] < z1[0]) return -1;
    if(z0[0] > z1[0]) return 1;

    
    
    if(zero[0] < one[0]) return -1;
    if(zero[0] > one[0]) return 1;

    if(zero[1] < one[1]) return -1;
    if(zero[1] > one[1]) return 1;

    if(zero[2] < one[2]) return -1;
    if(zero[2] > one[2]) return 1;
    
    return 0;
}