#include<conio.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<string.h>
#include<dos.h>


/******************************* set the following three  parameters ********/
/******************************* output  is the file solmaxg.txt ********/


/******* vpin is the rank of graphs generated ***********/
#define     vpin      9
//pow2=2^vpin
#define     pow2      512
//clexpect=3.2^([vpin/2]-2)+[vpin/2]-vpin
#define     clexpect   7




/****** mu is the eigenvalue, here it is the eigenvalue 0 *********/
#define     mu        0
#define     fwenable  1
#define     mxcnt     500
char        swg[100];
int         mtrxin[vpin+1][vpin+1][vpin+1];
int         mtrx[vpin+1][vpin+1];
int         mtrxasl[100][100];
int         mtrxcon[vpin+1][vpin+1];
int         mtrxinasl[vpin+1][vpin+1];
int         mtrxvec[vpin*vpin];
FILE        *fpg;
int         shfti[vpin*vpin];
int         shftj[vpin*vpin];
int         statvec[vpin*vpin];
int         last,indx,nosoltri,deti,cntasl,nomaxg,nosoltrifull,aut,clexpect1;
int         vecs[pow2][vpin+1];
int         plc[vpin+1];
int         savvecscomp[vpin+1];
int         savvecsgph[pow2+1];
int         noposlvl[pow2+1];
int         vallvl[pow2+1];
int         poslvl[mxcnt][mxcnt];
int         mtrxgph[mxcnt][mxcnt];
int         mtrxiso[100][100][100];
int         degiso[100][100];


/**************generates (0,1) vectors of length vpin  ***********/
void genvecs()
{
   int   i,j;

   plc[0]=1;
   for(i=1;i<=vpin;i++)
       plc[i]=plc[i-1]<<1;
   for(i=0;i<pow2;i++)
   {
       for(j=0;j<vpin;j++)
       if(i & plc[j])
          vecs[i][j+1]=1;
   }
}


void initial()
{
   int   i,j,d1=1;

   statvec[d1]=1;
   for(j=1;j<=vpin;j++)
   {
     for(i=1;i<j;i++)
     {
          shfti[d1]=i;
          shftj[d1]=j;
          d1++;
     }
     if(d1)
          statvec[d1-1]=1;
   }
   last=d1-1;
   mtrxvec[0]=1;
   indx=1;
   mtrxvec[indx]=2;
}


/****** tests if a generated graph of order vpin is in canonical form *********/
int  cantest(int np)
{
     int  i,j,e,curiso;
     int  prmiso[vpin+1];
     int  cprmiso[vpin+1];

     aut=0;
     curiso=1;
     prmiso[1]=0;
     for(i=1;i<=np;i++)
          cprmiso[i]=0;
     while(1)
     {
          curiso++;
          do
          {
                if(!(--curiso))
                        return(1);
                cprmiso[prmiso[curiso]]=0;
                for(i=prmiso[curiso]+1;i<=np;i++)
                if(!cprmiso[i])
                        break;
          }while(i>np);
          prmiso[curiso]=i;
          cprmiso[i]=1;
          for(j=1;j<curiso;j++)
          {
                if((e=mtrx[prmiso[j]][i]-mtrx[j][curiso])>0)
                        return(0);
                if(e<0)
                        break;
          }
          if(j==curiso)
          {
                if((++curiso)>np)
                {
                        curiso--;
                        aut++;
                }
                else
                        prmiso[curiso]=0;
          }
     }
}


void frd()
{
     int  i,j,sum;

     for(i=1;i<=vpin;i++)
     for(j=1;j<=vpin;j++)
          mtrxin[vpin][i][j]=mtrx[i][j];
     for(i=1;i<=vpin;i++)
          mtrxin[vpin][i][i]=-mu;
     for(i=1;i<=vpin;i++)
     for(j=1;j<=vpin;j++)
          mtrxinasl[i][j]=mtrxin[vpin][i][j];
     if(mu==0 || mu==-1)
     for(j=1;j<=vpin;j++)
     {
          sum=0;
          for(i=1;i<=vpin;i++)
          if(mtrxinasl[i][j])
              sum+=plc[i-1];
          savvecscomp[j]=sum;
     }
}


/***** tests if a generated graph is full rank ************/
int det(int lvl)
{
   int   i,j,k,l,ddd;

   if(lvl==2)
        ddd=mtrxin[lvl][1][1]*mtrxin[lvl][2][2]-mtrxin[lvl][1][2]*mtrxin[lvl][2][1];
   else
   {
        ddd=0;
        for(i=1;i<=lvl;i++)
        if(mtrxin[lvl][1][i])
        {
             for(j=2;j<=lvl;j++)
             {
                l=0;
                for(k=1;k<=lvl;k++)
                if(k!=i)
                {
                   l++;
                   mtrxin[lvl-1][j-1][l]=mtrxin[lvl][j][k];
                }
             }
             if(i%2)
                ddd=ddd+mtrxin[lvl][1][i]*det(lvl-1);
             else
                ddd=ddd-mtrxin[lvl][1][i]*det(lvl-1);
        }
   }
   return(ddd);
}


void frdspcl(int ii, int jj)
{
     int  i,j,l,k;

     k=0;
     for(i=1;i<=vpin;i++)
     if(i!=ii)
     {
        k++;
        l=0;
        for(j=1;j<=vpin;j++)
        if(j!=jj)
        {
             l++;
             mtrxin[vpin-1][k][l]=mtrxinasl[i][j];
        }
     }
}


/******* finds the inverse of a full rank graph of order vpin *******/
void  conversem()
{
   int i,j,x;

   for(i=1;i<=vpin;i++)
   for(j=1;j<=vpin;j++)
   {
          frdspcl(i,j);
          x=det(vpin-1);
          if((i+j)%2)
             mtrxcon[i][j]=-x;
          else
             mtrxcon[i][j]=x;
   }
}


void test()
{
   int  i,j,l,sum;

   for(i=1;i<=vpin;i++)
   for(j=1;j<=vpin;j++)
   {
       sum=0;
       for(l=1;l<=vpin;l++)
           sum+=(mtrxinasl[i][l]*mtrxcon[l][j]);
       if(i!=j && sum)
       {
           printf("\n  ERRRRRRRRRRRRRRR");
           getchar();
       }
       if(i==j && sum!=deti)
       {
           printf("\n  ERRRRRRRRRRRRRRR");
           getchar();
       }
   }
}


/******* finds the verices of compatibility graph ********/
void findvertcs()
{
   int i,j,cnt=0,sum,l,f;

   for(i=1;i<pow2;i++)
   {
       for(j=1;j<=vpin;j++)
       if((mu==-1 || mu==0) && i==savvecscomp[j])
          break;
       if(j>vpin)
       {
          for(l=1;l<=vpin;l++)
          {
            for(f=l+1;f<=vpin;f++)
            if(vecs[i][l] && vecs[i][f] && mtrx[l][f]==1)
                 break;
            if(f<=vpin)
                 break;
          }
          if(l>vpin)
          {
           sum=0;
           for(l=1;l<=vpin;l++)
           for(f=1;f<=vpin;f++)
           if(vecs[i][l] && vecs[i][f])
               sum+=mtrxcon[l][f];
           if(sum==((-mu)*deti))
           {
               cnt++;
               savvecsgph[cnt]=i;
           }
          }
       }
   }
   if(cnt>=mxcnt)
   {
        printf("\n ERRRRRRRRRRRRRRRRR mxcnt ..........");
        getchar();
        exit(1);
   }

   for(i=1;i<=cnt;i++)
   for(j=i+1;j<=cnt;j++)
   {
          sum=0;
          for(l=1;l<=vpin;l++)
          for(f=1;f<=vpin;f++)
          if(vecs[savvecsgph[i]][l] && vecs[savvecsgph[j]][f])
               sum+=mtrxcon[l][f];
          if(sum==deti  ||  sum==0)
          {
               mtrxgph[i][j]=abs(sum)+1;
               mtrxgph[j][i]=abs(sum)+1;
          }
          else
          {
               mtrxgph[i][j]=0;
               mtrxgph[j][i]=0;
          }
   }
   cntasl=cnt;
}


void makeg()
{
   int i,j,f,n=clexpect1,x;


   for(i=1;i<=vpin;i++)
   for(j=i+1;j<=vpin;j++)
   {
       mtrxasl[i+clexpect1][j+clexpect1]=mtrx[i][j];
       mtrxasl[j+clexpect1][i+clexpect1]=mtrx[i][j];
   }
   for(j=1;j<=n;j++)
   {
       x=savvecsgph[vallvl[j]];
       for(f=1;f<=vpin;f++)
       {
          mtrxasl[f+n][j]=vecs[x][f];
          mtrxasl[j][f+n]=vecs[x][f];
       }
   }
   for(i=1;i<=n;i++)
   for(j=i+1;j<=n;j++)
   {
       if(mtrxgph[vallvl[i]][vallvl[j]]==1)
       {
           mtrxasl[i][j]=0;
           mtrxasl[j][i]=0;
       }
       else
       {
           mtrxasl[i][j]=1;
           mtrxasl[j][i]=1;
       }
   }
}


/****** tests if a graph is bipartite *******/
int tstbipartite()
{
    int i,j,k,tab=0,n=clexpect1+vpin;
    int statver[100];

    for(i=1;i<=n;i++)
        statver[i]=0;
    while(1)
    {
      if(tab==0)
      {
        for(i=1;i<=n;i++)
        if(statver[i]==0)
        {
           statver[i]=1;
           break;
        }
        if(i>n)
           return(1);
      }
      for(i=1;i<=n;i++)
      if(statver[i]==0)
      {
        for(j=1;j<=n;j++)
        if(mtrxasl[i][j]==1 && statver[j]!=0)
           break;
        if(j<=n)
        {
           for(k=j+1;k<=n;k++)
           if(mtrxasl[i][k]==1 && statver[k] && statver[k]!=statver[j])
                break;
           if(k<=n)
                return(0);
           else
                statver[i]=-statver[j];
           tab=1;
           break;
        }
      }
      if(i>n)
          tab=0;
    }
}


/****** tests if a graph is triangle-free *******/
int trifree()
{

   int i,j,k,n=clexpect1+vpin;

   for(j=1;j<=n;j++)
   for(i=1;i<j;i++)
   if(mtrxasl[i][j]==1)
   for(k=1;k<i;k++)
   if(mtrxasl[k][j]==1  && mtrxasl[i][k]==1)
        return(0);
   return(1);
}


/******** writes solutions in a file  *******/
wtg()
{
    int i,j,n=clexpect1+vpin;

    fpg=fopen(swg,"a");
    fprintf(fpg,"\n\n%d ",nomaxg);
    for(i=1;i<=n;i++)
    {
      fprintf(fpg,"\n");
      for(j=1;j<=n;j++)
         fprintf(fpg," %d",mtrxasl[i][j]);
    }
    fclose(fpg);
}


/*********************** checks  graphs for isomorphism ***********************/
int  isofind()
{
    int  n=clexpect1+vpin;
    int  ii,i,j,e,curiso;
    int  prmiso[100];
    int  cprmiso[100];
    int  deg[100];

    for(i=1;i<=n;i++)
        deg[i]=0;
    for(i=1;i<=n;i++)
    for(j=1;j<=n;j++)
    if(mtrxasl[i][j])
        deg[i]++;
    for(ii=0;ii<nomaxg;ii++)
    {
       curiso=1;
       prmiso[1]=0;
       for(i=1;i<=n;i++)
          cprmiso[i]=0;
       while(1)
       {
          curiso++;
          do
          {
                if(!(--curiso))
                        goto l1;
                cprmiso[prmiso[curiso]]=0;
                for(i=prmiso[curiso]+1;i<=n;i++)
                if(!cprmiso[i] && deg[curiso]==degiso[ii][i])
                        break;
          }while(i>n);
          prmiso[curiso]=i;
          cprmiso[i]=1;
          for(j=1;j<curiso;j++)
          if(mtrxiso[ii][prmiso[j]][i]!=mtrxasl[j][curiso])
                break;
          if(j==curiso)
          {
                if((++curiso)>n)
                        return(0);
                else
                        prmiso[curiso]=0;
          }
       }
       l1: ;
    }
    for(i=1;i<=n;i++)
    for(j=1;j<=n;j++)
       mtrxiso[nomaxg][i][j]=mtrxasl[i][j];
    for(i=1;i<=n;i++)
       degiso[nomaxg][i]=deg[i];
    return(1);
}


int clique(int lvl,int val)
{
   int  i,j,k,x;

   vallvl[lvl]=poslvl[lvl][val];
   x=vallvl[lvl];

   for(i=1;i<lvl;i++)
   if(mtrxgph[vallvl[i]][x]>1)
   for(k=1;k<=vpin;k++)
   if(vecs[savvecsgph[vallvl[i]]][k]==1 && vecs[savvecsgph[x]][k]==1)
       return(1);

   for(i=1;i<lvl;i++)
   if(mtrxgph[vallvl[i]][x]>1)
   for(k=1;k<i;k++)
   if(mtrxgph[vallvl[k]][vallvl[i]]>1 && mtrxgph[vallvl[k]][x]>1)
       return(1);

   j=0;
   for(i=val+1;i<noposlvl[lvl];i++)
   if(mtrxgph[poslvl[lvl][i]][x])
       poslvl[lvl+1][j++]=poslvl[lvl][i];

   if(!j)
   {
       if(lvl==clexpect)
       {
           clexpect1=lvl;
           makeg();
           if(trifree()  && !tstbipartite())
           {
               printf("\n \n  TRIFREE..............................\n\n");
               if(isofind())
               {
                  wtg();
                  nomaxg++;
               }
           }
       }
       if(lvl>clexpect)
       {
           clexpect1=lvl;
           makeg();
           if(trifree()  && !tstbipartite())
           {
                printf("\n  ERRRRRROOORRRR!!! The theorem does not hold.........");
                getchar();
           }
       }
       return(1);
   }
   else
   {
       noposlvl[lvl+1]=j;
       for(i=0;i<noposlvl[lvl+1];i++)
          clique(lvl+1,i);
   }
   return(1);
}


/******************** finds  cliques in compatibility graph ***********/
int cliquem()
{
    int i,j;

    if(cntasl>=clexpect)
    {
            j=0;
            for(i=1;i<=cntasl;i++)
                poslvl[1][j++]=i;
            if(j)
            {
                noposlvl[1]=j;
                for(i=0;i<noposlvl[1];i++)
                    clique(1,i);
            }
    }
    return(1);
}


int  forward()
{
   int  i,j,ii,jj,k;

   mtrxvec[indx]--;
   ii=shfti[indx];
   jj=shftj[indx];
   mtrx[ii][jj]=mtrxvec[indx];
   mtrx[jj][ii]=mtrxvec[indx];
   if(!statvec[indx]  && mtrxvec[indx]==1 && mtrx[ii][jj-1]==0)
   {
       for(i=1;i<ii;i++)
       if(mtrx[i][jj]!=mtrx[i][jj-1])
         break;
       if(i==ii)
         return(2);
   }

   if(statvec[indx]  && jj<=(vpin))
   {
       for(i=1;i<jj;i++)
       if(mtrx[i][jj]==1)
       for(k=1;k<i;k++)
       if(mtrx[k][jj]==1  && mtrx[i][k]==1)
          return(2);
       if(!cantest(jj))
          return(2);
   }
   if(indx==last)
   {
       nosoltri++;
       if(nosoltri%1000==0)
         printf("\n  nosoltri=%d ",nosoltri);
       frd();
       deti=det(vpin);
       if(deti)
       {
           nosoltrifull++;
           if(fwenable)
           {
              conversem();
              test();
              findvertcs();
              cliquem();
           }
       }
       return(2);
   }
   indx++;
   mtrxvec[indx]=2;
   return(1);
}


int  back()
{
   while(mtrxvec[indx]==0)
        indx--;
   if(indx==0)
        return(0);
   return(1);
}


/************ a backtracking algorithm to generate triangle-free full rank graphs of order vpin ********/
void bk()
{
     int flg=1;
     do
     {
        switch(flg)
        {
            case 1:  flg=forward();
                     break;
            case 2:  flg=back();
                     break;
        }
     }while(flg);
}


main()
{
   genvecs();
   initial();
   sprintf(swg,"solmaxg.txt");
   fpg=fopen(swg,"w");
   fclose(fpg);
   bk();
   printf("\n\n\n ENNNDDDDDDDDD...\n\n  no  of  triangle-free graphs on %d vertices=%d \n\n  no  of full rank triangle-free graphs on %d vertices=%d \n\n  no  of maximum order non-bipartite triangle-free graphs of rank %d=%d\n\n",vpin,nosoltri,vpin,nosoltrifull,vpin,nomaxg);
   getchar();
}