#include <stdio.h>
#define S 20000 // how many binary ops to remember for isomorphism check
#define P 30  // how many transformations to remember for isomorphism check
#define Nlim 5 // 4, 5
#define outputBO 0 // 0, 1
#define remem isA // 0, isA, 1

#if outputBO
#define printLastC() {printf("Isomorphism class %d: T=%d, A=%d, C=%d\n",\
  I, C, a_class, c_class ); for(j=0;j<N;++j) { for(i=0;i<C;++i) { \
  for(k=0;k<N;++k) printf(" %d",f[h-C+i+1][j][k]); printf("  "); } \
  printf("\n"); } printf("\n"); }
#else
#define printLastC() {printf("Isomorphism class %d: T=%d, A=%d, C=%d\n",\
  I, C, a_class, c_class ); }
#endif

// associso.c, Christian van den Bosch, December 2002
// generate closed binary operations (BOs) on sets of size N.
// check each BO for associativity and commutativity.
// check pairs of binary operations for isomorphism by comparing
// them to each other's images under bijective transformation,
// and hence group them into isomorphism classes.

int main()
{
  int bo[4][4]; // current BO
  int f[S][4][4]; // list of BOs
  int lA[S], lC[S]; // assoc/commut flags for BOs.
  int p[4], q[4]; // transformation bijection and inverse
  int t[P][4];  // transformation bijection list
  int N; // size of set to investigate
  for(N=2; N<Nlim; N++)
  {
    long int a=0, b=0, c=0; // count of ABOs, bijections, CBOs
    long int d=0, o=0; // count of CABOs, all BOs
    int h, i, j, k, l; // counters
    int x,y,z; // three operands for associativity checking (base-N counter)
    int isA, done, isC; // booleans
    int I=0, C=0; // number of isomorphism classes, number in class
    int m=0; // how many BOs to consider for isomorphism
    int isOk; // boolean
    int a_class=0, c_class=0;
    // start binary operations from 0
    for( i=0; i<N; ++i )
      for( j=0; j<N; ++j )
        bo[i][j]=0;
    // now, start looking at all possible closed BOs
    printf( "Investigating closed binary operations on set of order %d...\n", N);
    i = j = 0;
    while(i<N && j<N) // done when counter overflows
    {
      // is this a CBO (commutative BO)?
      isC = 1; // assume so for now
      for( i=0; i<N; ++i )
        for( j=i+1; j<N; ++j )
          if( bo[i][j] != bo[j][i] ) // if i*j != j*i
            isC = 0; // not a CBO
      // is this an ABO (associative BO)?
      isA=1; // assume so for now
      done=x=y=z=0; // start operand counter at 000
      while(isA && !done) // check until not an ABO, or counter overflows
      { // not an ABO if x*(y*z) != (x*y)*z
        isA &= ( bo[x][bo[y][z]] == bo[bo[x][y]][z]);
        // increment xyz representing our three operands
        ++x; // increment, and
        if(!(x%=N)) // if zero, then carry
        {
          ++y; // increment, and
          if(!(y%=N)) // if zero, then carry
          {
            ++z; // increment, and
            if(!(z%=N)) // if zero, then overflow
              done=1; // meaning we've checked all possible operands
          }
        }
      }
      if( remem ) // save BOs for isomorphism testing with other BOs
      {
        for( i=0; i<N; ++i )
          for( j=0; j<N; ++j )
            f[m][i][j] = bo[i][j];
        lC[m]=isC; // remember this BO's commutativity
        lA[m]=isA; // remember this BO's associativity
        ++m; // remember how many BO's we're checking for isomorphism
      }
      ++o; // increase total BO count
      c += isC; // if commutative, increase CBO count
      a += isA; // if associative, increase ABO count
      d += ( isC & isA ); // if both, increase CABO count
      // increment N*N digit base-N counter representing BO
      for(i=0; i<N; ++i)
      {
        for(j=0; j<N; ++j)
        {
          ++bo[i][j]; // increment digit
          if( bo[i][j]%=N ) // if zero, loop and carry.
          { // if non-zero, don't loop and carry.
            goto allCarried; // ugly code alert!
          } // but them's the breaks
        } // when all you really want
      } // is a double break.
allCarried:
    } // while( i<N && j<N ) - fall out of loop when counter overflows
    // report what we've found.
    printf( "Counts for order %d: T=%ld, ", N, o );
    printf( "A=%ld, C=%ld, CA=%ld\n\n", a, c, d );
    // now generate bijective functions for isomorphism checking.
    // we use an N-digit base-N counter p to generate functions.
    // each bijective function is stored for later use.
    // we begin with the non-bijective function where p[i]=0 for all i.
    for( i=0; i<N; ++i )
      p[i]=0;
    b=i=0; // we have no bijections yet.
    while( i<N ) // loop until the counter overflows
    {
      // function is a bijection if it's invertible.
      // write rogue values to the inverse function
      // to ease detection of non invertible functions
      for( i=0; i<N; ++i )
         q[i]=-1;
      // write the inverse of the function
      for( i=0; i<N; ++i )
         q[p[i]] = i;
      // check that every part of the inverse was written precisely once.
      isOk = 1; // innocent until proven guilty.
      for( i=0; i<N; ++i )
      {
        if( q[p[i]] != i ) // if anything was overwritten or unwritten,
        {
          isOk = 0;      // p is not invertible, hence not a bijection.
          break;
        }
      }
      if( isOk ) // p invertible, so is a bijection.
      {
        // we don't want the identity bijection.
        isOk = 0; // assume the worst
        for( i=0; i<N; ++i )
          if( p[i] != i) // if any p[i] != i
            isOk = 1; // it's not the identity
        if( isOk ) // if it's not the identity transformation
        {
          // we want it, so save it
          for( i=0; i<N; ++i )
            t[b][i]=p[i];
          // and its inverse, in case we want it later.
          ++b; // increment bijection count
        }
      }
      // increment counter representing function,
      // falling out of while loop when it overflows.
      for( i=0; i<N; ++i )
      {
        ++p[i]; // increment digit
        if( p[i]%=N ) // if zero, loop and carry.
          break; // if non-zero, don't loop and carry
      }
    } // while( i<N ) - all functions checked when counter overflows.
    // there are much more efficient ways to generate transformations,
    // but this one is nice and general, and allows some code reuse.

    // now compare BOs to each others' transformations under
    // bijection, in order to find isomorphism classes
    isOk = 0; // we haven't found any isomorphisms yet
    for( h=0; h<m; ++h ) // step thru BO list
    {
      // if we haven't just found an isomorphism and printed it out,
      // print the current BO. (if we had, this would be it.)
      if( !isOk )
      { // print the current BO
        C=1;
        a_class = lA[h];
        c_class = lC[h];
	++I;
      }
      // look for BO f[i] ( o ) isomorphic to BO f[h] ( * )
      for( i=h+1; i<m; ++i )
      {
        // check if * and o (f[h] and f[i]) are isomorphic under alpha.
        for( l=0; l<b; ++l ) // step through list of bijections (alpha)
        {
          // now check isomorphism of * and o under alpha for all j, k
          isOk = 1; // assume isomorphic under this alpha for now
          for( j=0; j<N; ++j ) // j is left operand to BO
          {
            for( k=0; k<N; ++k ) // k is right operand to BO
            { // compare (j*k)alpha to (j)alpha o (k)alpha
              if( t[l][f[h][j][k]] != f[i][t[l][j]][t[l][k]] ) // if unequal,
              {
                isOk = 0; // * and o are not isomorphic under this alpha
                goto nextAlpha; // so try the next alpha
              }
            }
          }
          if( isOk ) // if * and o are isomorphic under alpha,
            break; // don't look for another alpha
nextAlpha:
        } // now isOk tells us whether * and o are isomorphic under alpha.
        if( isOk )
        {
          a_class += (l = lA[i]);
          lA[i] = lA[h+1];
          lA[h+1] = l;
          c_class += (l = lC[i]);
          lC[i] = lC[h+1];
          lC[h+1] = l;
          // swap f[i] with f[h+1], element by element
          for( j=0; j<N; ++j )
          {
            for( k=0; k<N; ++k )
            {
              l = f[i][j][k];
              f[i][j][k] = f[h+1][j][k];
              f[h+1][j][k] = l;
            }
          }
          ++C;
          break; // stop looking for . - we found one!
        }
      } // for( i ) ( o )
      // any o found to match * becomes the new *
      // any other new * will be of a different isomorphism class
      if( !isOk )
        printLastC();
    } // for( h ) ( * )
    if( isOk )
      printLastC();
  }
}