#include<stdlib.h>
#include<stdio.h>
#include<math.h>

/*This algorithm aims to compute the homogeneous boolean rank and it is based on this rank's connection to the Johnson graph. The relation between these two concepts is explained in "On Ranks of Regular Polygons".
Given n and k, it determines if it is possible to construct a homogeneous boolean factorization of size k of the boolean matrix of n-gons.
If so, such factorization is returned; otherwise, two zero matrices are shown.
The least integer for which a homogeneous boolean factorization of the boolean matrix of n-gons exists is the homogeneous boolean rank of n-gons.

/*Sub-functions*/

/*Given a vector of integers, a, of size n, this function checks if all its entries are distinct*/
int unique(int a[], int n)
{
  int i, k;

  for (i = 0; i < n; i++) {
    for (k = i+1; k < n; ++k) {
        if(a[i]==a[k]){
          return 0;
        }
    }
  }
  return 1;
}

/* Function that checks if a given integer, value, belongs to a specific integer vector (sorted in ascending order), x, of size n*/
int mymember(int x[], int value, int n)
{

    int t=-1, first = 0, last=n-1, middle = (first+last)/2;
    while (first <= last)
    {
        if (x[middle] < value)
        {
           first = middle+1;
        }
        else
        {
            if (x[middle] == value)
            {
                t=middle;
                break;
            }
            else
            {
                last=middle-1;
            }
        }
        middle = (first + last)/2;
    }
    return t;

}

/* Function that returns the value of "n choose k" */
int C(int n, int k)
{
    if ( k==0 || k==n )
    {
        return 1;
    }
    else
    {
        return C(n-1,k-1) + C(n-1,k);
    }
}

/*Function that determines if the cardinality of the symmetric difference of two zero-one vectors is two*/
int xor2(int v[], int u[], int n)
{
    int sum=0, i=0;
    while(sum<=2 && i<n)
    {
        sum=sum+v[i]+u[i]-2*v[i]*u[i];
        i=i+1;
    }
    return (sum==2);
}

/* Function that determines if two zero-one vectors, v and u, both of size n, are equal */
int arequal(int v[], int u[], int n)
{
    int i=0;
    while((v[i]==u[i])&&(i<n))
        i=i+1;
    if(i==n)
        return 1;
    else
        return 0;
}


/*Subprogram used to generate the indicator vectors of all combinations of M elements drawn without replacement
  from a set of N elements (as implemented in http://www.netlib.org/toms-2014-06-10/382)

Reference:
Phillip J Chase, `Algorithm 382: Combinations of M out of N Objects [G6]',
Communications of the Association for Computing Machinery 13:6:368 (1970).*/

int twiddle(x, y, z, p)
int *x, *y, *z, *p;
  {
  register int i, j, k;
  j = 1;
  while(p[j] <= 0)
    j++;
  if(p[j-1] == 0)
    {
    for(i = j-1; i != 1; i--)
      p[i] = -1;
    p[j] = 0;
    *x = *z = 0;
    p[1] = 1;
    *y = j-1;
    }
  else
    {
    if(j > 1)
      p[j-1] = 0;
    do
      j++;
    while(p[j] > 0);
    k = j-1;
    i = j;
    while(p[i] == 0)
      p[i++] = -1;
    if(p[i] == -1)
      {
      p[i] = p[k];
      *z = p[k]-1;
      *x = i-1;
      *y = k-1;
      p[k] = -1;
      }
    else
      {
      if(i == p[0])
	return(1);
      else
	{
	p[j] = p[i];
	*z = p[i]-1;
	p[i] = 0;
	*x = j-1;
	*y = i-1;
	}
      }
    }
  return(0);
  }

  void inittwiddle(m, n, p)
int m, n, *p;
  {
  int i;
  p[0] = n+1;
  for(i = 1; i != n-m+1; i++)
    p[i] = 0;
  while(i != n+1)
    {
    p[i] = i+m-n;
    i++;
    }
  p[n+1] = -2;
  if(m == 0)
    p[1] = 1;
  }

/*Main function*/
void main()
{
  int n, k, m;
  printf("Value of n:");
  scanf("%d", &n);
  printf("Value of k:");
  scanf("%d", &k);
  m=k/2;
  int i, x, y, z, p[k+2], b[k];

  /* Set of possible path vertices, represented by a matrix, V. Each line of V represents the indicator vector of a certain vertex*/
  inittwiddle(m, k, p);
  int C1= C(k,m);
  int V[C1][k];


  for(i = 0; i != k-m; i++)
    {
    b[i] = V[0][i]=0;
    }

  while(i != k)
    {
    V[0][i]=1;
    b[i++]=1;
    }

  int j=0;

  while(!twiddle(&x, &y, &z, p))
    {
    j++;
    b[x] = 1;
    b[y] = 0;
    for(i = 0; i != k; i++)
      {
      if(b[i])
        V[j][i]=1;
      else
        V[j][i]=0;
      }
    }

    /* Set of possible path edge colors, represented by a matrix, Col. Each line of Col is the indicator vector of an edge color*/
    inittwiddle(k-(m+1), k, p);
    int C2=C(k,k-(m+1));
    int Col[C2][k];



  for(i = 0; i != (m+1); i++)
    {
    b[i] = Col[0][i]=0;
    }

  while(i != k)
    {
    Col[0][i]=1;
    b[i++]=1;
    }

  j=0;

  while(!twiddle(&x, &y, &z, p))
    {
    j++;
    b[x] = 1;
    b[y] = 0;
    for(i = 0; i != k; i++)
      {
      if(b[i])
        Col[j][i]=1;
      else
        Col[j][i]=0;
      }
    }

    /*Construction of matrix Adj, whose line i contains the vertices adjacent to vertex V[i]*/
    int count;
    int Adj[C1][m*(k-m)];

    for(i=0;i<C1;i++)
    {
    count=-1;
    for(j=0;j<C1;j++)
    {
        if(xor2(V[i],V[j],k))
         {
            count=count+1;
            Adj[i][count]=j;
         }
    }
    }

    /*Construction of two matrices:
    - CAdj, whose entry ij is the color of the edge that joins V[i] to Adj[i][j];
    - Incv, whose line i contains all vertices that are incident to an edge of color Col[i]*/
    int n_adj=m*(k-m);
    int CAdj[C1][n_adj];
    int countcol[C2];
    for(i=0;i<C2;i++)
    {
        countcol[i]=0;
    }
    int Incv[C2][m+1];
    for(i=0;i<C2;i++)
    {for(j=0;j<m+1;j++)
        {Incv[i][j]=-1;
        }
    }

    int symdif[k];

    int l;

    for(i=0;i<C1;i++)
    {
        for(j=0;j<n_adj;j++)
        {
        count=0;
        for(l=0;l<k;l++)
        {symdif[l]=1-((V[i][l]+V[Adj[i][j]][l])>0);
        }

        while(!(arequal(symdif,Col[count],k)))
        {
         count=count+1;
        }

        CAdj[i][j]=count;
        if (!(Incv[count][countcol[count]]==i))
        {
        if (Incv[count][countcol[count]]>-1)
         {countcol[count]=countcol[count]+1;}
         Incv[count][countcol[count]]=i;}

        }

    }

    /*Sorting of CAdj lines (ascending order)*/
    int temp1,temp2;
    for(l=0;l<C1;l++)
    {
    for(i = 0; i < n_adj; i++)
    {
    for(j = i; j < n_adj; j++)
      {
      if(CAdj[l][j] < CAdj[l][i])
       {
        temp1 = CAdj[l][i];
        temp2 = Adj[l][i];
        CAdj[l][i] = CAdj[l][j];
        Adj[l][i] = Adj[l][j];
        CAdj[l][j] = temp1;
        Adj[l][j] = temp2;
       }
       }
    }
    }

    /*Path/Cycle construction

    P[i], CP[i] and sizep represent the path ith vertex, the color of the path ith edge and the path size, respectively.

    The vector ind is used to do the depth-first search, through the adjacent vertices of each vertex.*/

    int P[n+1], CP[n], ind[C1], sizep=1;
    for(i=0; i<(n+1); i++)
    {
        P[i]=0;
    }
    for(i=0;i<C1;i++)
    {
        ind[i]=0;
    }

    /*Symmetry classes

    If CD^T is a homogeneous boolean factorization of the boolean matrix of an n-gon, the same goes for (C*P)*(D*P)^T, where P is a permutation matrix.
    In that sense, we consider that two paths/cycles lie on the same symmetry class if the corresponding homogeneous boolean factorizations are, up to permutations, equal.
    This means that we can improve the algorithm by skipping paths that lie on symmetry classes that have already been analyzed.
    To do that, we assign an integer to each column of the first matrix of the factorization (each column can be seen as a binary number). That integer is just the conversion of the binary number to decimal system.
    Thus, for each factorization, we get a set of k integers.
    In the end, two paths are on the same symmetry class if their sets of k integers are the same.


    Matrix Sym stores the symmetry classes.
    size_sym stores the number of symmetry classes for the different path sizes*/

    int Sym[n][n_adj][k], size_sym[n], aux[k], aux_ord[k], Sym_ord[n][n_adj][k];

    size_sym[0]=1;
    for(j=0;j<k;j++)
    {
        Sym[0][0][j]=V[0][j];
    }
    for(i=1;i<n;i++)
    {
        size_sym[i]=0;
    }


    /*Boolean variables that turn on/off symmetry analysis*/
    int bool_sym=1, stop=0,

    int c=1, q, d;


    /*Loop used to construct the path*/
    while(c && sizep<=n)
    {
        /*Possibility of closing the path*/
        if(sizep==n)
        {
            d=0;
            for(i=0;i<n_adj;i++)
            {
                if (P[0]=Adj[P[n-1]][i])
                {
                    d=1;
                    l=i;
                }
            }
            i=2;
            while(d && i<sizep)
            {
                /*Restrictions on the edge colors and on the vertices*/
                d =!((mymember(CAdj[P[i-1]],CAdj[P[n-1]][l],n_adj)+1) || (mymember(Incv[CP[i-1]],P[0],m+1)+1));
                i=i+1;
            }
            /*If the restrictions are satisfied, the path is augmented (and we get a cycle)*/
            if(d)
            {
                sizep=sizep+1;
                P[n]=P[0];
                CP[n-1]=CAdj[P[n-1]][l];
            }
            /* Otherwise, it is decreased*/
            else
            {
                ind[P[n-1]]=0;
                sizep=sizep-1;
                ind[P[sizep-1]]=ind[P[sizep-1]]+1;
            }
        }
        else
        {
            d=1;
            j=ind[P[sizep-1]];
            while(d && j<n_adj)
            {
                i=1;
                while(d && i<sizep)
                {
                    /*Restrictions on the edge colors and on the vertices*/
                    d=!((mymember(CAdj[P[i-1]],CAdj[P[sizep-1]][j],n_adj)+1) || (mymember(Incv[CP[i-1]], Adj[P[sizep-1]][j],m+1)+1));
                    i=i+1;
                }
                /*Restrictions satisfied*/
                if(d)
                {
                    /*Symmetry Analysis on*/
                    if(bool_sym)
                    {
                        /*Set of k integers associated with the path*/
                        for(i=0;i<k;i++)
                        {
                            aux[i]=Sym[sizep-1][size_sym[sizep-1]-1][i]+(((int)pow(2,sizep))*V[Adj[P[sizep-1]][j]][i]);
                        }
                        for(i=0;i<k;i++)
                        {
                            aux_ord[i]=aux[i];
                        }
                        for(l=0;l<k;l++)
                        {
                            for(i = l; i < k; i++)
                            {
                                 if(aux_ord[i] < aux_ord[l])
                                {
                                    temp1 = aux_ord[l];
                                    aux_ord[l] = aux_ord[i];
                                    aux_ord[i]=temp1;
                                }
                            }
                        }
                        /*Does the current path  belong to any of the symmetry classes already analyzed?*/
                        i=0;
                        while(i<size_sym[sizep])
                        {
                            l=0;
                            while((l<k) && (aux_ord[l]==Sym_ord[sizep][i][l]))
                            {
                                l=l+1;
                            }
                            if(l==k)
                            {
                                break;
                            }
                            else
                            {
                                i=i+1;
                            }
                        }

                        /*If it doesn't, a new symmetry class is added*/
                        if(i==size_sym[sizep])
                        {
                            size_sym[sizep+1]=0;
                            size_sym[sizep]=size_sym[sizep]+1;
                            for(i=0;i<k;i++)
                            {
                                Sym[sizep][size_sym[sizep]-1][i]=aux[i];
                                Sym_ord[sizep][size_sym[sizep]-1][i]=aux_ord[i];
                            }
                            /*Path augmentation*/
                            ind[P[sizep-1]]=j;
                            CP[sizep-1]=CAdj[P[sizep-1]][j];
                            P[sizep]=Adj[P[sizep-1]][j];
                            sizep=sizep+1;
                            d=0;
                            if(unique(Sym[sizep-2][size_sym[sizep-2]-1],k))
                            {
                                bool_sym=0;
                                stop=sizep-2;
                            }
                        }
                        /*If it does, there is no need to analyze it. The search goes on.*/
                        else
                        {
                            j=j+1;
                        }
                    }
                    else
                    {
                        /*Path augmentation*/
                        ind[P[sizep-1]]=j;
                        CP[sizep-1]=CAdj[P[sizep-1]][j];
                        P[sizep]=Adj[P[sizep-1]][j];
                        d=0;
                        sizep=sizep+1;
                    }

                }
                else
                {
                    /*The search proceeds */
                    j=j+1;
                    d=1;
                }

            }
            /*None of the adjacent vertices meets the restrictions. The path is decreased*/
            if(j>=n_adj)
            {
                ind[P[sizep-1]]=0;
                sizep=sizep-1;
                c=(sizep>=1);
                /*Symmetry Analysis is reactivated*/
                if(sizep==stop)
                {
                    bool_sym=1;
                }
                if(c)
                {
                    ind[P[sizep-1]]=ind[P[sizep-1]]+1;
                }
            }
        }
    }

/*Zero matrices*/

int C[n][k], D[n][k];

for(i=0;i<n;i++)
{
    for(j=0;j<k;j++)
        {
            C[i][j]=0;
            D[i][j]=0;
        }
}

/*If the path has been built, the corresponding homogeneous boolean factorization is constructed*/
if(sizep>2)
{
    for(i=0;i<n;i++)
    {
        for(j=0;j<k;j++)
        {
            C[i][j]=V[P[i]][j];
            D[i][j]=Col[CP[i]][j];
        }
    }

}

/*OUTPUT*/

FILE *fp;

char buf[0x100];
snprintf(buf,sizeof(buf),"OUTPUTn%dk%d.txt",n,k);

fp=fopen(buf,"w");

fprintf(fp,"OUTPUT \nn=%d \n",n);
fprintf(fp,"k=%d \n \n",k);

fprintf(fp,"C=[");
for(i = 0; i<n-1; i++)
{
    for(j = 0; j<k; j++)
    {
        fprintf(fp,"%d ", C[i][j]);
    }
    fprintf(fp,"\n");
}
for(j = 0; j<k; j++)
    {
        fprintf(fp,"%d ", C[n-1][j]);
    }
    fprintf(fp,"]\n \n");

fprintf(fp,"D=[");
for(i = 0; i<n-1; i++)
{
    for(j = 0; j<k; j++)
    {
        fprintf(fp,"%d ", D[i][j]);
    }
    fprintf(fp,"\n");
}
for(j = 0; j<k; j++)
    {
        fprintf(fp,"%d ", D[n-1][j]);
    }
    fprintf(fp,"]\n");
}