/*                        isom.c                     11-11-89 

     This contains all the procedures pertaining to the
     elimination of isomorphic copies from the output of
     matrix-generating programs such as MaGIC.c.
      
     It is assumed that the external variables siz, neg[], 
     ord[][], des and C[][] are available and that this 
     module is called, presumably by the accept() function,
     every time a new matrix is found.  If there are any 
     isomorphic copies, the integer isoms is incremented.
      
     There are two parts to all this.  First come some
     procedures to find the acceptable permutations on the
     current siz, neg[], ord[][] and des.  Then there are
     the routines to be executed on discovery of a new 
     matrix.                                               */


/***********************************************************
***********************************************************/

perm_initial()              /*** Set up perm and istak ****/
{ register int i;

  perm = (PRM*) malloc(sizeof(PRM));
  perm->h[0] = -1;
  perm->pup = 0;
  
  for (i=0; i<ISOMAX; i++)
  { **(istak[i].ic) = '\0';
    istak[i].left = istak[i].right = istak[i].parent = 0;
  }
}


/***********************************************************
***********************************************************/

setperm()     /*** Find the automorphisms which do not 
                   disturb ord, N or des.              ***/
{ register int i,j;
  int p_poss[SZ][SZ];
  int p_used[SZ], sh[SZ], single[SZ];
  PRM *prmptr;
  
  for ( prmptr = perm; prmptr->h[0] >= 0; prmptr = prmptr->pup )
  prmptr->h[0] = -1;
  FORALL(i)  FORALL(j)
  {
     p_poss[i][j] = (lower_than(i)==lower_than(j)) &&
                    (higher_than(i)==higher_than(j));
     if ( theJob->f_n && 
          (( N[i]==i && N[j]!=j ) || ( N[j]==j && N[i]!=i )) )
     p_poss[i][j] = 0;
  }
  FORALL(i)
  { 
     if ( theJob->f_t )
     {
        if ( i != des ) 
        p_poss[des][i] = p_poss[i][des] = 0;
        if ( theJob->f_n && i != N[des] ) 
        p_poss[N[des]][i] = p_poss[i][N[des]] = 0;
     }
     else if ( true[i] ) 
     FORALL(j) if ( !true[j] ) 
     p_poss[i][j] = p_poss[j][i] = 0;
     p_used[i] = 1;  sh[i] = i;
     single[i] = 1;
     FORALL(j) if ( i != j && p_poss[i][j] ) single[i] = 0;
  }
  goto CHANGE;               /** !! **/
  
PTEST:  FORALL(i)  FORALL(j)
        if ( ord[i][j] != ord[sh[i]][sh[j]] ) goto CHANGE;
        newperm(sh);
CHANGE: FORALL(i) if ( !(( theJob->f_n && i < N[i] ) || single[i] ))
        { 
           p_used[sh[i]] = 0;
           if ( theJob->f_n ) p_used[sh[N[i]]] = 0;
           for ( j = sh[i]-1; j > 0; j-- )
           if ( p_poss[i][j] && !p_used[j] )
           { 
              sh[i] = j;  if ( theJob->f_n ) sh[N[i]] = N[j];
              p_used[j] = 1;
              if ( theJob->f_n ) p_used[N[j]] = 1;
              for ( i--; i >= 0 && !(theJob->f_n && i < N[i]); i-- ) 
              if ( !single[i] )
              { 
                 for (j = siz; j; j-- )
                 if ( p_poss[i][j] && !p_used[j] )
                 { 
                    sh[i] = j;  if ( theJob->f_n ) sh[N[i]] = N[j];
                    p_used[j] = 1;
                    if ( theJob->f_n) p_used[N[j]] = 1;
                    break;
                 }
              }
              goto PTEST;
           }
        }
}


          /***********************/

lower_than(x)  int x;     /*** How many elements below x ***/
{ register int i;
  int j=0;

  for ( i=0; i < x; i++ ) j += ord[i][x];
  return(j);
}


          /***********************/

higher_than(x)  int x;    /*** How many elements above x ***/
{ register int i;
  int j=0;

  for ( i = siz; i > x; i-- ) j += ord[x][i];
  return(j);
}



          /**********************/

newperm(vec)  int *vec;           /*** Add a permutation ***/
{ int i;
  PRM *p;
  
  for ( p = perm; p->pup && p->h[0] >= 0; p = p->pup) ;
  if ( !p->pup )
  {
     p->pup = (PRM*) malloc(sizeof(PRM));
     p->pup->pup = 0;
     p->pup->h[0] = -1;
  }
  FORALL(i) p->h[i] = *(vec+i);
}



/***********************************************************
***********************************************************/

isomorphic(ptr)   ISM *ptr;   /*** Matrix was already in the
                              "isom" tree below ptr. If so,
                              snip it out; if not add its 
                              copies.                   ***/
{ register int i,j;
  
  FORALL(i) FORALL(j)
  { 
     if ( C[i][j] < ptr->ic[i][j] )
     { 
        if ( ptr->left ) 
        return(isomorphic(ptr->left));
        add_isoms();
        return(0);
     }
     if ( C[i][j] > ptr->ic[i][j] )
     { 
        if ( ptr->right ) 
        return(isomorphic(ptr->right));
        add_isoms();
        return(0);
     }
  }
  snip(ptr);
  isoms++;
  return(1);
}



          /******************************/

snip(p)  ISM *p;         /*** Remove *p from "isom" tree ***/
{ ISM *q;
  int i,j;

  if ( p == istak )
  { 
     if ( !p->left && !p->right ) 
     **(p->ic) = 0;
     else
     { 
        if ( p->right ) 
        for ( q = p->right; q->left; q = q->left ) ;
        else
        for ( q = p->left; q->right; q = q->right ) ;
        FORALL(i) FORALL(j)
        p->ic[i][j] = q->ic[i][j];
        snip(q);
     }
  }
  else
  { 
     if ( !p->left ) subst(p,p->right);
     else if ( !p->right ) subst(p,p->left);
     else
     { 
        for ( q = p->right; q->left; q = q->left ) ;
        FORALL(i) FORALL(j)
        p->ic[i][j] = q->ic[i][j];
        subst(q,q->right);
     }
  }
}


          /******************************/

subst(p1,p2)  ISM *p1,*p2;  /*** Move p2 into place of p1.
                                 Mark p1 as defunct      ***/
{
  if ( p2 ) p2->parent = p1->parent;
  if ( p1->parent->left == p1 ) p1->parent->left = p2;
  else p1->parent->right = p2;
  **(p1->ic) = 0;
}



          /******************************/

add_isoms()        /*** Put nontrivial isomorphic copies 
                        of C[][] into the "isom" tree.   ***/
{ int i,j,same;
  char ac[SZ][SZ];
  PRM *aptr;
  
  for ( aptr = perm; aptr->h[0] >= 0; aptr = aptr->pup )
  { 
     FORALL(i) FORALL(j) 
     ac[aptr->h[i]][aptr->h[j]] = aptr->h[C[i][j]];
     same = 1;
     FORALL(i) FORALL(j)
     if ( ac[i][j] != C[i][j] )
     { 
        same = 0;
        i = j = siz;
     }
     if ( !same ) 
     { 
        if ( **((*(istak)).ic) ) add_this(ac,istak);
        else
        { 
           FORALL(i) FORALL(j)
           (*(istak)).ic[i][j] = ac[i][j];
        }
     }
  }
}



          /******************************/

add_this(mat,i_tree)  char mat[SZ][SZ]; ISM *i_tree;
                        /*** If matrix "mat" not in tree
                             below "i_tree", put it in.  ***/
{ register int i,j;
  ISM *tack_on();

  FORALL(i) FORALL(j)
  if ( mat[i][j] < i_tree->ic[i][j] )
  { 
     if ( i_tree->left ) return(add_this(mat,i_tree->left));
     i_tree->left = tack_on(i_tree,mat);
     return(1);
  }
  else if ( mat[i][j] > i_tree->ic[i][j] )
  { 
     if ( i_tree->right ) return(add_this(mat,i_tree->right));
     i_tree->right = tack_on(i_tree,mat);
     return(1);
  }
  return(0);
}



          /*********************************/

ISM *tack_on(p,mat)  ISM *p; char mat[SZ][SZ];
                          /*** New ism with parent p
                               and matrix mat        ***/
{ int i,j,k;

  for ( i = 0; i < ISOMAX; i++ )
  if ( !**(istak[i].ic) )
  { 
     istak[i].left = istak[i].right = 0;
     istak[i].parent = p;
     FORALL(j) FORALL(k)
     istak[i].ic[j][k] = mat[j][k];
     return(istak+i);
  }
  skipout("Isomorphism stack overflow");
}