/* Analyzer for automorphisms, vertex equivalence classes,
   Also includes code to produce input to AT&T graph formatters */
/* Bisection width calculation is based on dynamic programming
   technique in J.A. Lukes, Combinatorial Solution to the
   Partitioning of General Graphs, IBM J. on R & D,
   Vol 19, March 1975.
   Coded recursively to keep storage management simple */
/* Bisection width calculation performed by METIS */
/* Graph formatter inputs are written to separate files */
/* Adjacency lists replace adjacency matrices */

#include <stdarg.h>             /* va_start, va_arg, va_end */
#include <stdlib.h>             /* EXIT_FAILURE, malloc, qsort */
#include <string.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <math.h>
#define MAXN 32675
#include "../NAUTY/nauty.h"
#include "../metis-3.0/Lib/defs.h"
#include "../metis-3.0/Lib/struct.h"
/* #include "../metis-3.0/Lib/macros.h" not needed, conflicts w/ nauty.h*/
#include "../metis-3.0/Lib/rename.h"
#include "../metis-3.0/Lib/proto.h"

/* FORMATTING is to indicate that an execution is to produce
   input to neato (undirected) or dot (directed).
   INCLUDELABELS indicates that labels are desired in nodes
   ANALYZING indicates that analysis information is to be written to
   stdout.
   APPROXBISECTION indicates the largest graph size (in vertices) for
   which Lukes's algorithm should be used.  If larger than this, METIS
   is used.
   EDGESGUESS is a guess at the number of edges per vertex.
   DISCARDSELFLOOPS indicates whether reflexive (i,i) edges are removed.
*/
#define FORMATTING TRUE
#define INCLUDELABELS FALSE
#define ANALYZING TRUE
#define APPROXBISECTION 40
#define EDGESGUESS 20
#define DISCARDSELFLOOPS FALSE

/* For NAUTY */
graph g[MAXN*MAXM];
nvector lab[MAXN],ptn[MAXN],orbits[MAXN];
static DEFAULTOPTIONS(options);
statsblk(stats);
setword workspace[50*MAXM];
int numWords,v;
set *gv;
/*****************/


FILE *scriptDescript;  /* script for running formatters */
float nautyCPU=0.0,metisCPU=0.0,memoryCPU=0.0,hashCPU=0.0,LukesCPU=0.0,
  permutationCPU=0.0,sortCPU;
int numVertices,numEdges,vertexCount;
struct edge {
  int tail,head;
};
typedef struct edge edgeType;
edgeType edgeTab[EDGESGUESS*MAXN];
int *tailHeader,*headTab,*headHeader,*tailTab;
char **vertexLabel;
int *hashLabel;
char graphName[60];

/* For Lukes's bisection width computation */
int W;
int **conn,**adjacency;
int *workPartition;
int partitionCount,*bestPartition,bestValue;

void memoryDie()
{
printf("Out of memory\n");
exit(0);
}

float CPUseconds()
{
struct rusage rusage;

getrusage(RUSAGE_SELF,&rusage);
return rusage.ru_utime.tv_sec+rusage.ru_utime.tv_usec/1000000.0;
}

int twoPower(k)
/* Computes powers of two as int */
int k;
{
return 1<<k;
}

double twoPowerBig(k)
/* Computes powers of two as double */
int k;
{
int i;
double power=1;

for (i=0;i<k;i++)
  power*=2.0;
return power;
}

void allocate()
/* Allocates tables that receive extended use */
{
int i,j;
float startCPU;

startCPU=CPUseconds();
numEdges=0;  /* Perhaps generators should set & then verify? */
vertexLabel=(char**) malloc(numVertices*sizeof(char*));
if (!vertexLabel)
  memoryDie();
hashLabel=(int*) malloc(2*numVertices*sizeof(int));
if (!hashLabel)
  memoryDie();
for (i=0;i<2*numVertices;i++)
  hashLabel[i]=(-1);
vertexCount=0;
if (numVertices<=APPROXBISECTION)
{
  conn=(int**) malloc(numVertices*sizeof(int*));
  adjacency=(int**) malloc(numVertices*sizeof(int*));
  if (!conn || !adjacency)
    memoryDie();
  for (i=0;i<numVertices;i++)
  {
    conn[i]=(int*) malloc((numVertices+1)*sizeof(int));
    adjacency[i]=(int*) malloc((numVertices+1)*sizeof(int));
    if (!conn[i] || !adjacency[i])
      memoryDie();
  }
  workPartition=(int*) malloc(numVertices*sizeof(int));
  if (!workPartition)
    memoryDie();
}
bestPartition=(int*) malloc(numVertices*sizeof(int));
if (!bestPartition)
  memoryDie();
memoryCPU+=CPUseconds()-startCPU;
}

void deallocate()
/* Deallocates memory claimed in allocate() */
{
int i;
float startCPU;

startCPU=CPUseconds();
for (i=0;i<vertexCount;i++)
  free(vertexLabel[i]);
free(vertexLabel);
free(hashLabel);
free(tailHeader);
free(headTab);
free(headHeader);
free(tailTab);
if (numVertices<=APPROXBISECTION)
{
  for (i=0;i<numVertices;i++)
  {
    free(conn[i]);
    free(adjacency[i]);
  }
  free(conn);
  free(adjacency);
  free(workPartition);
}
free(bestPartition);
memoryCPU+=CPUseconds()-startCPU;
}

int label2number(char *string)
/* Finds label in linear probe hash table.  If not present, then insert */
{
int i,j;
float startCPU;

startCPU=CPUseconds();
i=0;
for (j=0;string[j]!='\0';j++)
  i+=abs(string[j]);
i=i%(2*numVertices);
while (hashLabel[i]!=(-1) && strcmp(string,vertexLabel[hashLabel[i]])!=0)
  i=(i+1)%(2*numVertices);
hashCPU+=CPUseconds()-startCPU;
if (hashLabel[i]!=(-1))
  return hashLabel[i];
/*
printf("New label: %s\n",string);
*/
hashLabel[i]=vertexCount;
if (vertexCount==numVertices)
{
  printf("Too many vertex labels\n");
  abort();
}
vertexLabel[vertexCount]=strdup(string);
vertexCount++;
return hashLabel[i];
}

void connectVertices(char* tail,char* head)
{
int i,j;

if (numEdges==EDGESGUESS*MAXN)
{
  printf("Too many edges\n");
  exit(0);
}
i=label2number(tail);
j=label2number(head);
if (i==j && DISCARDSELFLOOPS)
{
  printf("Warning - self-loop removed\n");
  return;
}
edgeTab[numEdges].tail=i;
edgeTab[numEdges].head=j;
numEdges++;
}

int tailThenHead(const void* xin, const void* yin)
{
int result;
edgeType *x,*y;

x=(edgeType*) xin;
y=(edgeType*) yin;
result=x->tail - y->tail;
if (result!=0)
  return result;
else
  return x->head - y->head;
}

int headThenTail(const void* xin, const void* yin)
{
int result;
edgeType *x,*y;

x=(edgeType*) xin;
y=(edgeType*) yin;
result=x->head - y->head;
if (result!=0)
  return result;
else
  return x->tail - y->tail;
}

void buildAdjacencies()
{
int duplicates=0,i,j,k;
float startCPU;

tailHeader=(int*) malloc((numVertices+1)*sizeof(int));
headTab=(int*) malloc(numEdges*sizeof(int));
headHeader=(int*) malloc((numVertices+1)*sizeof(int));
tailTab=(int*) malloc(numEdges*sizeof(int));
if (!tailHeader || !headTab || !headHeader || !tailTab)
  memoryDie();
startCPU=CPUseconds();
qsort(edgeTab,numEdges,sizeof(edgeType),tailThenHead);
sortCPU+=CPUseconds()-startCPU;
j=k=0;
for (i=0;i<numVertices;i++)
{
  tailHeader[i]=k;
  while (j<numEdges && edgeTab[j].tail==i)
  {
    if (j>0 && edgeTab[j].tail==edgeTab[j-1].tail &&
      edgeTab[j].head==edgeTab[j-1].head)
    {
      duplicates++;
      j++;
      continue;
    }
    headTab[k++]=edgeTab[j++].head;
  }
}
tailHeader[numVertices]=k;
if (duplicates>0)
  printf("%d duplicates found & removed\n",duplicates);

startCPU=CPUseconds();
qsort(edgeTab,numEdges,sizeof(edgeType),headThenTail);
sortCPU+=CPUseconds()-startCPU;
j=k=0;
for (i=0;i<numVertices;i++)
{
  headHeader[i]=k;
  while (j<numEdges && edgeTab[j].head==i)
  {
    if (j>0 && edgeTab[j].tail==edgeTab[j-1].tail &&
      edgeTab[j].head==edgeTab[j-1].head)
    {
      j++;
      continue;
    }
    tailTab[k++]=edgeTab[j++].tail;
  }
}
headHeader[numVertices]=k;

numEdges-=duplicates;
}

void connectedSet()
{
/* See p. 172 of Lukes' paper for definition */
/* Warshall's algorithm has been adapted */
int i,j,k,**work;  /* work is the number of vertices
                      on the path */

work=(int**) malloc(numVertices*sizeof(int*));
if (!work)
  memoryDie();
for (i=0;i<numVertices;i++)
{
  work[i]=(int*) malloc(numVertices*sizeof(int));
  if (!work[i])
    memoryDie();
}

for (i=0;i<numVertices;i++)
  for (k=0;k<numVertices;k++)
    work[i][k]=99999;
for (i=0;i<numVertices;i++)
  for (k=tailHeader[i];k<tailHeader[i+1];k++)
    work[i][headTab[k]]=2;

for (j=numVertices-1;j>=0;j--)
{
  for (i=0;i<numVertices;i++)
    if (work[i][j]!=99999)
      for (k=0;k<numVertices;k++)
        if (work[j][k]!=99999)
          if (work[i][j]+work[j][k]-1<work[i][k])
            work[i][k]=work[i][j]+work[j][k]-1;
/*
  printf("conn[%d]=",j);
*/
  k=0;
  for (i=0;i<j;i++)
    if (work[i][j]<=W)
    {
      conn[j][k++]=i;
/*
      printf("%d ",i);
*/
    }
  conn[j][k]=(-1);
/*
  printf("\n");
*/
}
for (i=0;i<numVertices;i++)
{
  k=0;
  for (j=tailHeader[i];j<tailHeader[i+1];j++)
    adjacency[i][k++]=headTab[j];
  adjacency[i][k]=99999;
}

for (i=0;i<numVertices;i++)
  free(work[i]);
free(work);
}

void LukesAlgorithm(nextVertex,inputValue,
  inputLeftWeight,inputRightWeight,remainingEdges)
int nextVertex,inputValue,
  inputLeftWeight,inputRightWeight,remainingEdges;
{
int i,workValue,workRemainingEdges;

if (nextVertex==numVertices)
{
  partitionCount++;
  if (inputValue>bestValue)
  {
/**/
    printf("bisection improved from %d to %d\n",bestValue,inputValue);
/**/
    bestValue=inputValue;
    for (i=0;i<numVertices;i++)
      bestPartition[i]=workPartition[i];
  }
  return;
}
if (inputValue+remainingEdges<=bestValue)
  return;  /* Cut-off search - this subtree will not contain global optimum */
if (inputLeftWeight<W)
{
  for (i=0;
       conn[nextVertex][i]!=(-1) && workPartition[conn[nextVertex][i]]==1;
       i++)
    ;
  if (conn[nextVertex][i]!=(-1))
  {
    workPartition[nextVertex]=0;
    workValue=inputValue;
    workRemainingEdges=remainingEdges;
    for (i=0;adjacency[nextVertex][i]<nextVertex;i++)
    {
      workRemainingEdges--;
      if (workPartition[adjacency[nextVertex][i]]==0)
        workValue++;
    }
    LukesAlgorithm(nextVertex+1,workValue,inputLeftWeight+1,inputRightWeight,
      workRemainingEdges);
  }
}
if (inputRightWeight==0)
{
  workPartition[nextVertex]=1;
  workRemainingEdges=remainingEdges;
  for (i=0;adjacency[nextVertex][i]<nextVertex;i++)
    workRemainingEdges--;
  LukesAlgorithm(nextVertex+1,inputValue,inputLeftWeight,inputRightWeight+1,
    workRemainingEdges);
  return;
}
if (inputRightWeight<W)
{
  for (i=0;
       conn[nextVertex][i]!=(-1) && workPartition[conn[nextVertex][i]]==0;
       i++)
    ;
  if (conn[nextVertex][i]!=(-1))
  {
    workPartition[nextVertex]=1;
    workValue=inputValue;
    workRemainingEdges=remainingEdges;
    for (i=0;adjacency[nextVertex][i]<nextVertex;i++)
    {
      workRemainingEdges--;
      if (workPartition[adjacency[nextVertex][i]]==1)
        workValue++;
    }
    LukesAlgorithm(nextVertex+1,workValue,inputLeftWeight,inputRightWeight+1,
      workRemainingEdges);
  }
}
}

void nautyVerify(str,computedAutos1,computedAutos2,computedEquivClasses,
  computedDiameter,computedIndegree,computedOutdegree,
  computedBisectionWidth)
char *str;
double computedAutos1;
int computedAutos2,computedEquivClasses,computedDiameter,
  computedIndegree,computedOutdegree,computedBisectionWidth;
{
int i,j,k;
char computedAutos1string[100],grpsize1string[100];
int currentAttribute=0,*nodeAttributes;
char *colorTable[7]={"red","blue","green","yellow","cyan","violet",
  "orange"};
int indegree,outdegree,diameter,bisectionWidth;
int *queue,tail,head,*dist,*mark;
int wgtflag=0,Cstyle=0,nparts=2,defaultOptions=0;
int partitionSize[2];
char fileName[40];
char neatoExtension[5]=".nto",dotExtension[5]=".dot";
FILE *fileDescript;
char *labelString,emptyLabel[10]="label=\"\",",defaultLabel[1]="";
char *nodeHeightString,smallHeight[11]="height=0.1",
  largeHeight[11]="height=0.5";
char *nodeWidthString,smallWidth[10]="width=0.1",largeWidth[10]="width=0.5";
float startCPU;

buildAdjacencies();

printf("%s ",str);
printf("%d vertices ",numVertices);
if (options.digraph)
  printf("directed\n");
else
  printf("undirected\n");
fflush(stdout);
if (vertexCount!=numVertices)
{
  if (numVertices==1)
    return;
  printf("%s: missing vertices!!!\n",str);
  return;
}
options.writemarkers=FALSE;
options.writeautoms=FALSE;
options.digraph=FALSE;
if (!DISCARDSELFLOOPS)
{
  /* Check for self-loops */
  for (i=0;i<numVertices;i++)
  {
    for (j=tailHeader[i];
         j<tailHeader[i+1] && headTab[j]<i;
         j++)
      ;
    if (j<tailHeader[i+1] && headTab[j]==i)
      options.digraph=TRUE;
  }
}
for (i=0;i<=numVertices;i++)
  if (headHeader[i]!=tailHeader[i])
  {
    options.digraph=TRUE;
    break;
  }
if (!options.digraph)
{
  for (i=0;i<numEdges;i++)
    if (headTab[i]!=tailTab[i])
    {
      options.digraph=TRUE;
      break;
    }
}

if (FORMATTING)
{
  for (i=0;i<39 && str[i]!=')';i++)
    if (str[i]==',' || str[i]=='(')
      fileName[i]='.';
    else
      fileName[i]=str[i];
  if (options.digraph)
    for (k=0;i+k<39 && k<4;k++)
      fileName[i+k]=dotExtension[k];
  else
    for (k=0;i+k<39 && k<4;k++)
      fileName[i+k]=neatoExtension[k];
  fileName[i+k]='\0';
  if ((fileDescript=fopen(fileName,"w"))==NULL)
  {
    printf("file %s open failed\n",fileName);
    return;
  }
  if (options.digraph)
    fprintf(scriptDescript,"dot -Tps %s -o /dev/null\n",fileName);
  else
    fprintf(scriptDescript,"neato -Gepsilon=.2 -v -Tps %s -o /dev/null\n",
      fileName);
}

numWords = (numVertices + WORDSIZE - 1) / WORDSIZE;
for (i=0;i<numVertices;i++)
{
  gv = GRAPHROW(g,i,numWords);
  EMPTYSET(gv,numWords);
  for (j=headHeader[i];j<headHeader[i+1];j++)
    ADDELEMENT(gv,tailTab[j]);
}

startCPU=CPUseconds();
nauty(g,lab,ptn,NILSET,orbits,&options,&stats,workspace,
  50*MAXM,numWords,numVertices,NILGRAPH);
nautyCPU+=CPUseconds()-startCPU;

/* bisection information is computed here so that it
   is available when neato is used */
if (options.digraph)
  bisectionWidth=0;
else
{
  if (numVertices<=APPROXBISECTION) /* Lukes code */
  {
    startCPU=CPUseconds();
    W=(numVertices+1)/2;
    connectedSet();
    workPartition[0]=0;
    partitionCount=0;
    bestValue=0;
/*  numEdges is now computed during tail-adjacency table generation */
    LukesAlgorithm(1,0,1,0,numEdges/2);
    LukesCPU+=CPUseconds()-startCPU;
    if (ANALYZING)
      printf("complete partitions reached %d\n",partitionCount);
    if (partitionCount==0)
      bisectionWidth=(-1);
    else
      bisectionWidth=numEdges/2-bestValue;
/*
    printf("best partition has value %d\n",bestValue);
    printf("left partition: ");
    for (i=0;i<numVertices;i++)
      if (bestPartition[i]==0)
        printf("%d ",i);
    printf("\nright partition: ");
    for (i=0;i<numVertices;i++)
      if (bestPartition[i]==1)
        printf("%d ",i);
    printf("\n");
    printf("Edges that are cut: %d\n",numEdges/2-bestValue);
*/
  }
  if (numVertices>APPROXBISECTION || bisectionWidth==(-1))
  /* METIS code */
  {
    startCPU=CPUseconds();
    METIS_PartGraphRecursive(&numVertices,
                             tailHeader,headTab,
                             NULL,NULL,
                             &wgtflag,
                             &Cstyle,
                             &nparts,
                             &defaultOptions,
                             &bisectionWidth,
                             bestPartition);
    metisCPU+=CPUseconds()-startCPU;
    partitionSize[0]=partitionSize[1]=0;
    for (i=0;i<numVertices;i++)
      partitionSize[bestPartition[i]]++;
    if (abs(partitionSize[0]-partitionSize[1])>1 && ANALYZING)
      printf("METIS returned unbalanced partitions %d %d!!!\n",
        partitionSize[0],partitionSize[1]);
    queue=(int*) malloc(numVertices*sizeof(int));
    mark=(int*) malloc(numVertices*sizeof(int));
    if (!queue || !mark)
      memoryDie();
    for (i=0;i<numVertices;i++)
      mark[i]=FALSE;
    tail=head=0;
    for (i=0;i<numVertices && bestPartition[i]!=0;i++)
      ;
    if (i<numVertices)
    {
      mark[i]=TRUE;
      queue[tail++]=i;
    }
    while (tail!=head)
    {
      i=queue[head++];
      for (j=tailHeader[i];j<tailHeader[i+1];j++)
        if (bestPartition[headTab[j]]==0 && !mark[headTab[j]])
        {
          mark[headTab[j]]=TRUE;
          queue[tail++]=headTab[j];
        }
    }
    for (i=0;i<numVertices && bestPartition[i]!=1;i++)
      ;
    if (i<numVertices)
    {
      mark[i]=TRUE;
      queue[tail++]=i;
    }
    while (tail!=head)
    {
      i=queue[head++];
      for (j=tailHeader[i];j<tailHeader[i+1];j++)
        if (bestPartition[headTab[j]]==1 && !mark[headTab[j]])
        {
          mark[headTab[j]]=TRUE;
          queue[tail++]=headTab[j];
        }
    }
    if (tail!=numVertices && ANALYZING)
      printf("METIS returned disconnected cluster!!!\n");
    free(queue);
    free(mark);
  }
}

if (FORMATTING)
{
  /* Code to produce dot (directed) or neato (undirected) input */
  nodeAttributes=(int*) malloc(numVertices*sizeof(int));
  if (!nodeAttributes)
    memoryDie();
  for (i=0;i<numVertices;i++)
    nodeAttributes[i]=(-1);
  if (options.digraph)
    fprintf(fileDescript,"di");
  fprintf(fileDescript,"graph G{\n");
  fprintf(fileDescript,"  font=helvetica;\n");
  fprintf(fileDescript,"  orientation=portrait;\n");
  fprintf(fileDescript,"  size=\"7.5,10\";\n");
  if (stats.numorbits>1)
  {
    /* Vertex display properties based on vecs */
    for (i=0;i<numVertices;i++)
    {
      if (nodeAttributes[orbits[i]]==(-1))
        nodeAttributes[orbits[i]]=currentAttribute++;
      nodeAttributes[i]=nodeAttributes[orbits[i]];
    }
  }
  else
  {
    /* Vertex display properties based on distance from vertex 0 */
    queue=(int*)malloc(numVertices*sizeof(int));
    if (!queue)
      memoryDie();
    tail=head=0;
    nodeAttributes[0]=0;
    queue[tail++]=0;
    while (tail!=head)
    {
      i=queue[head++];
      for (j=tailHeader[i];j<tailHeader[i+1];j++)
        if (nodeAttributes[headTab[j]]==(-1))
        {
          nodeAttributes[headTab[j]]=nodeAttributes[i]+1;
          queue[tail++]=headTab[j];
        }
    }
    free(queue);
  }
  for (i=0;i<numVertices;i++)
  {
    if (INCLUDELABELS)
    {
      labelString=defaultLabel;
      nodeWidthString=largeWidth;
      nodeHeightString=largeHeight;
    }
    else
    {
      labelString=emptyLabel;
      nodeWidthString=smallWidth;
      nodeHeightString=smallHeight;
    }
/*  Shapes only, no color
    fprintf(fileDescript,
      "\"%s\" [%sshape=polygon,sides=%d,peripheries=%d,font=helvetica,\n",
      vertexLabel[i],labelString,nodeAttributes[i]%5+4,
      nodeAttributes[i]/5+1);
    fprintf(fileDescript,"  %s,%s];\n",nodeWidthString,nodeHeightString);
*/
/*  Colors only, no shapes
    fprintf(fileDescript,
      "\"%s\" [%sfont=helvetica,\n",
      vertexLabel[i],labelString);
    fprintf(fileDescript,
      "  color=%s,style=filled,\n",colorTable[nodeAttributes[i]%7]);
    fprintf(fileDescript,"  %s,%s];\n",nodeWidthString,nodeHeightString);
*/
/*  Shapes and colors - corrected to avoid excessive peripheries
    Colors change most rapidly, then # of sides, then # of peripheries
*/
    fprintf(fileDescript,
      "\"%s\" [%sshape=polygon,sides=%d,peripheries=%d,font=helvetica,\n",
      vertexLabel[i],labelString,(nodeAttributes[i]/7)%5+4,
      nodeAttributes[i]/35+1);
    fprintf(fileDescript,
      "  color=%s,style=filled,\n",colorTable[nodeAttributes[i]%7]);
    fprintf(fileDescript,"  %s,%s];\n",nodeWidthString,nodeHeightString);
  }
  free(nodeAttributes);
  if (options.digraph)
  {
    for (i=0;i<numVertices;i++)
      for (j=tailHeader[i];j<tailHeader[i+1];j++)
        fprintf(fileDescript,
          "\"%s\"->\"%s\";\n",vertexLabel[i],vertexLabel[headTab[j]]);
  }
  else
  {
  /* Color the partitions as blue and green, but the bisecting
     edges are colored red */
    for (i=0;i<numVertices;i++)
    {
      for (j=tailHeader[i];
           j<tailHeader[i+1] && headTab[j]<=i;
           j++)
        ;
      for (;j<tailHeader[i+1];j++)
      {
        fprintf(fileDescript,
          "\"%s\"--\"%s\" ",vertexLabel[i],vertexLabel[headTab[j]]);
        if (bestPartition[i]!=bestPartition[headTab[j]])
          fprintf(fileDescript,"[color=red];\n");
        else if (bestPartition[i]==0)
          fprintf(fileDescript,"[color=blue];\n");
        else /* bestPartition[i]==1 */
          fprintf(fileDescript,"[color=green];\n");
      }
    }
  }
  fprintf(fileDescript,"}\n");
  fclose(fileDescript);
}
if (!ANALYZING)
  return;

indegree=0;
outdegree=0;
for (i=0;i<numVertices;i++)
{
  k=tailHeader[i+1]-tailHeader[i];
  if (k>outdegree)
    outdegree=k;
  k=headHeader[i+1]-headHeader[i];
  if (k>indegree)
    indegree=k;
}

/* Run breadth-first search once for each vec */
queue=(int*) malloc(numVertices*sizeof(int));
dist=(int*) malloc(numVertices*sizeof(int));
if (!queue || !dist)
  memoryDie();
diameter=0;
for (i=0;i<numVertices;i++)
  if (orbits[i]==i) /* Found a vec */
  {
    for (j=0;j<numVertices;j++)
      dist[j]=99999;
    tail=head=0;
    dist[i]=0;
    queue[tail++]=i;
    while (tail!=head)
    {
      j=queue[head++];
      if (dist[j]>diameter)
        diameter=dist[j];
      for (k=tailHeader[j];k<tailHeader[j+1];k++)
        if (dist[j]+1<dist[headTab[k]])
        {
          dist[headTab[k]]=dist[j]+1;
          queue[tail++]=headTab[k];
        }
    }
  }
free(queue);
free(dist);

/*Normalize automorphism counts*/
while (computedAutos2<stats.grpsize2)
{
  computedAutos2++;
  computedAutos1/=10;
}
if (computedAutos1==stats.grpsize1
    && computedAutos2==stats.grpsize2
    && computedEquivClasses==stats.numorbits
    && diameter==computedDiameter
    && indegree==computedIndegree
    && outdegree==computedOutdegree
    && bisectionWidth==computedBisectionWidth)
{
  printf("verified\n");
  return;
}

if (computedAutos1==0.0)
  printf("autos (nty) %fE%d\n",stats.grpsize1,stats.grpsize2);
else if (computedAutos1==stats.grpsize1
     && computedAutos2==stats.grpsize2)
  printf("autos verified\n");
else
{
  sprintf(computedAutos1string,"%f",computedAutos1);
  sprintf(grpsize1string,"%f",stats.grpsize1);
/*
  if (strcmp(computedAutos1string,grpsize1string)==0 &&
     computedAutos2==stats.grpsize2)
*/
  if (fabs(computedAutos1-stats.grpsize1)/computedAutos1<0.000001)
    printf("autos verified within accuracy\n");
  else
    printf("autos mismatch: computed=%fE%d nauty=%fE%d\n",
      computedAutos1,computedAutos2,stats.grpsize1,stats.grpsize2);
}

if (computedEquivClasses==0)
  printf("v.e.c. (nty) %d\n",stats.numorbits);
else if (computedEquivClasses==stats.numorbits)
  printf("v.e.c. verified\n");
else
  printf("v.e.c. mismatch: computed=%d nauty=%d\n",
    computedEquivClasses,stats.numorbits);

if (computedDiameter==0)
  printf("diameter (alg) %d\n",diameter);
else if (diameter==computedDiameter)
  printf("diameter verified\n");
else
  printf("diameter mismatch: computed=%d algorithm=%d\n",
    computedDiameter,diameter);

if (computedIndegree==0)
  printf("in-degree (alg) %d\n",indegree);
else if (indegree==computedIndegree)
  printf("in-degree verified\n");
else
  printf("in-degree mismatch: computed=%d algorithm=%d\n",
    computedIndegree,indegree);

if (computedOutdegree==0)
  printf("out-degree (alg) %d\n",outdegree);
else if (outdegree==computedOutdegree)
  printf("out-degree verified\n");
else
  printf("out-degree mismatch: computed=%d algorithm=%d\n",
    computedOutdegree,outdegree);

if (computedBisectionWidth==0)
  printf("bisection width (alg) %d\n",bisectionWidth);
else if (bisectionWidth==computedBisectionWidth)
  printf("bisection width verified\n");
else
  printf("bisection width mismatch: computed=%d algorithm=%d\n",
    computedBisectionWidth,bisectionWidth);
}

void readInput()
{
int i,j,numEdges;
char tail[30],head[30];

scanf("%d %d",&numVertices,&numEdges);
allocate();
for (i=0;i<numEdges;i++)
{
  scanf("%s %s",tail,head);
  connectVertices(tail,head);
}
nautyVerify("file input",0.0,0,0,0,0,0,0);
deallocate();
}

int powerCount(int bits,int radixLeaf)
{
/* Computes the power of two needed to get
   number of automorphisms for a binary heap tree
*/

int i;
for (i=0;(radixLeaf>>i)&1==1;i++)
  ;
if (i==0)
  i=1;
return radixLeaf-bits+i+(twoPower(bits)-1-radixLeaf)/2;

}

int vecCount(int bits,int radixLeaf)
{
/* Computes the number of vertex equivalence
   classes for a binary heap tree
*/
int work;
int i;

for (i=0;(radixLeaf>>i)&1==1;i++)
  ;
work=i+1+bits*(bits+1)/2-i*(i+1)/2;
for (;i<bits;i++)
  if ((radixLeaf>>i)&1==1)
    work++;
return work;
}

void tree(n)
int n;
{
int i,j,bits,leafNumber;
double autos;
int vec,indegree,outdegree,diameter;
char string1[20],string2[20];

sprintf(graphName,"tree(%d)",n);
numVertices=n;
allocate();
for (i=2;i<=n;i++)
{
  j=i/2;
  sprintf(string1,"%d",i);
  sprintf(string2,"%d",j);
  connectVertices(string1,string2);
  connectVertices(string2,string1);
}
bits=0;
leafNumber=n-1;
while (leafNumber>=twoPower(bits))
{
  leafNumber-=twoPower(bits);
  bits++;
}
if (n==1)
{
  autos=1.0;
  vec=1;
  indegree=0;
  outdegree=0;
}
else if (n==2)
{
  autos=2.0;
  vec=1;
  indegree=1;
  outdegree=1;
}
else if (n==3 || n==4)
{
  autos=2.0;
  vec=2;
  indegree=2;
  outdegree=2;
}
else
{
  autos=twoPowerBig(powerCount(bits,leafNumber));
  vec=vecCount(bits,leafNumber);
  indegree=3;
  outdegree=3;
}
diameter=0;
for (i=2;i<=n;i*=2)
  diameter++;
i/=2;
diameter*=2;
if (n<i+i/2)
  diameter--;
nautyVerify(graphName,autos,0,vec,diameter,indegree,outdegree,0);
deallocate();
}

void linearCBT(numTrees,treeHeight)
int numTrees,treeHeight;
{
int i,j,nodesInOneTree,leavesInOneTree,internalNodesInOneTree;
double autos1;
int autos2,vecs,diameter,indegree,outdegree,bisectionWidth;
char string1[25],string2[25];

sprintf(graphName,"linearCBT(%d,%d)",numTrees,treeHeight);
nodesInOneTree=twoPower(treeHeight+1)-1;
leavesInOneTree=twoPower(treeHeight);
internalNodesInOneTree=nodesInOneTree-leavesInOneTree;
numVertices=numTrees*nodesInOneTree;
allocate();
for (i=1;i<=numTrees;i++)
{
  for (j=2;j<=internalNodesInOneTree;j++)
  {
    sprintf(string1,"tree%dinternal%d",i,j);
    sprintf(string2,"tree%dinternal%d",i,j/2);
    connectVertices(string1,string2);
    connectVertices(string2,string1);
  }
  if (leavesInOneTree>1)
    for (j=1;j<=leavesInOneTree;j++)
    {
      sprintf(string1,"tree%dleaf%d",i,j);
      sprintf(string2,"tree%dinternal%d",i,(internalNodesInOneTree+j)/2);
      connectVertices(string1,string2);
      connectVertices(string2,string1);
    }
}
for (i=1;i<numTrees;i++)
  for (j=1;j<=leavesInOneTree;j++)
  {
    sprintf(string1,"tree%dleaf%d",i,j);
    sprintf(string2,"tree%dleaf%d",i+1,j);
    connectVertices(string1,string2);
    connectVertices(string2,string1);
  }
diameter=2*treeHeight+numTrees-1;
if (numTrees==1)
{
  autos1=twoPowerBig(twoPower(treeHeight)-1);
  autos2=0;
  vecs=treeHeight+1;
  if (treeHeight==0)
    indegree=outdegree=0;
  else if (treeHeight==1)
    indegree=outdegree=2;
  else
    indegree=outdegree=3;
}
else if (numTrees==2)
  if (treeHeight==0)
  {
    autos1=2.0;
    autos2=0;
    vecs=1;
    indegree=outdegree=1;
  }
  else if (treeHeight==1)
  {
    autos1=12.0;
    autos2=0;
    vecs=1;
    indegree=outdegree=2;
  }
  else
  {
    autos1=twoPowerBig(twoPower(treeHeight));
    autos2=0;
    vecs=treeHeight+1;
    indegree=outdegree=3;
  }
else
{
  autos1=twoPowerBig(twoPower(treeHeight));
  autos2=0;
  vecs=(treeHeight+1)*((numTrees+1)/2);
  if (treeHeight==0)
    indegree=outdegree=2;
  else
    indegree=outdegree=3;
}
if (numTrees==0)
  bisectionWidth=1;
else if (numTrees<leavesInOneTree)
  bisectionWidth=numTrees;
else
  bisectionWidth=leavesInOneTree;
nautyVerify(graphName,autos1,autos2,vecs,diameter,indegree,
  outdegree,bisectionWidth);
deallocate();
}

void meshCBT(rowTreeHeight,colTreeHeight)
int rowTreeHeight,colTreeHeight;
{
int i,j,nodesInOneRowTree,leavesInOneRowTree,nodesInAllRowTrees;
int internalNodesInOneRowTree;
int nodesInOneColTree,leavesInOneColTree,internalNodesInOneColTree;
int ijLeafNumber,ijColTreeParent;
double autos;
int vecs,diameter,indegree,outdegree,bisectionWidth;
char string1[20],string2[20];

if (rowTreeHeight>colTreeHeight)
  abort();
sprintf(graphName,"meshCBT(%d,%d)",rowTreeHeight,colTreeHeight);
nodesInOneRowTree=(twoPower(rowTreeHeight+1))-1;
leavesInOneRowTree=twoPower(rowTreeHeight);
internalNodesInOneRowTree=nodesInOneRowTree-leavesInOneRowTree;
nodesInOneColTree=(twoPower(colTreeHeight+1))-1;
leavesInOneColTree=twoPower(colTreeHeight);
internalNodesInOneColTree=nodesInOneColTree-leavesInOneColTree;
nodesInAllRowTrees=leavesInOneColTree*nodesInOneRowTree;
numVertices=nodesInAllRowTrees
  +leavesInOneRowTree*internalNodesInOneColTree;
allocate();
/* Row tree edges */
for (i=1;i<=leavesInOneColTree;i++)
{
  for (j=2;j<=internalNodesInOneRowTree;j++)
  {
    sprintf(string1,"row%dinternal%d",i,j);
    sprintf(string2,"row%dinternal%d",i,j/2);
    connectVertices(string1,string2);
    connectVertices(string2,string1);
  }
  if (leavesInOneRowTree>1)
    for (j=1;j<=leavesInOneRowTree;j++)
    {
      sprintf(string1,"row%dcol%d",i,j);
      sprintf(string2,"row%dinternal%d",i,(internalNodesInOneRowTree
        +j)/2);
      connectVertices(string1,string2);
      connectVertices(string2,string1);
    }
}
/* Column tree edges */
for (i=1;i<=leavesInOneRowTree;i++)
{
  for (j=2;j<=internalNodesInOneColTree;j++)
  {
    sprintf(string1,"col%dinternal%d",i,j);
    sprintf(string2,"col%dinternal%d",i,j/2);
    connectVertices(string1,string2);
    connectVertices(string2,string1);
  }
  if (leavesInOneColTree>1)
    for (j=1;j<=leavesInOneColTree;j++)
    {
      sprintf(string1,"row%dcol%d",j,i);
      sprintf(string2,"col%dinternal%d",i,(internalNodesInOneColTree
        +j)/2);
      connectVertices(string1,string2);
      connectVertices(string2,string1);
    }
}

if (rowTreeHeight==0)
{
  autos=twoPowerBig(twoPower(colTreeHeight)-1);
  vecs=colTreeHeight+1;
  if (colTreeHeight<2)
    indegree=outdegree=2;
  else
    indegree=outdegree=3;
}
else if (rowTreeHeight==colTreeHeight)
  if (rowTreeHeight==1)
  {
    autos=16.0;
    vecs=1;
    indegree=outdegree=2;
  }
  else
  {
    autos=twoPowerBig(twoPower(rowTreeHeight+1)-1);
    vecs=rowTreeHeight+1;
    indegree=outdegree=3;
  }
else
{
  autos=twoPowerBig(twoPower(rowTreeHeight)+twoPower(colTreeHeight)-2);
  vecs=rowTreeHeight+colTreeHeight+1;
  indegree=outdegree=3;
}
diameter=2*(rowTreeHeight+colTreeHeight);
bisectionWidth=twoPower(rowTreeHeight);
nautyVerify(graphName,autos,0,vecs,diameter,indegree,outdegree,
  bisectionWidth);
deallocate();
}

void CBTmesh(k)
int k;
{
int level,i,j,rows,columns;
int treeNode;
char string1[20],string2[20];

sprintf(graphName,"CBTmesh(%d)",k);
if (k==0)
  return;
numVertices=(2*k+1)*twoPower(k)*twoPower(k);
allocate();
/* Build mesh for each tree node */
rows=columns=twoPower(k);
for (level=0;level<2*k;level++)
{
  for (treeNode=twoPower(level);treeNode<twoPower(level+1);treeNode++)
    for (i=0;i<rows;i++)
      for (j=0;j<columns;j++)
      {
        sprintf(string1,"n%dr%dc%d",treeNode,i,j);
        if (i>0)
        {
          sprintf(string2,"n%dr%dc%d",treeNode,i-1,j);
          connectVertices(string1,string2);
        }
        if (i<rows-1)
        {
          sprintf(string2,"n%dr%dc%d",treeNode,i+1,j);
          connectVertices(string1,string2);
        }
        if (j>0)
        {
          sprintf(string2,"n%dr%dc%d",treeNode,i,j-1);
          connectVertices(string1,string2);
        }
        if (j<columns-1)
        {
          sprintf(string2,"n%dr%dc%d",treeNode,i,j+1);
          connectVertices(string1,string2);
        }
      }
  if (rows==columns)
    rows=rows/2;
  else
    columns=rows;
}
/* Connect for tree edges between levels */
rows=columns=twoPower(k);
for (level=0;level<2*k;level++)
{
  for (treeNode=twoPower(level);treeNode<twoPower(level+1);treeNode++)
  {
    /* Connect West to North for left side */
    for (i=0;i<rows;i++)
    {
      sprintf(string1,"n%dr%dc%d",treeNode,i,0);
      sprintf(string2,"n%dr%dc%d",2*treeNode,0,i);
      connectVertices(string1,string2);
      connectVertices(string2,string1);
    }
    /* Connect East to North for right side */
    for (i=0;i<rows;i++)
    {
      sprintf(string1,"n%dr%dc%d",treeNode,i,columns-1);
      sprintf(string2,"n%dr%dc%d",2*treeNode+1,0,rows-1-i);
      connectVertices(string1,string2);
      connectVertices(string2,string1);
    }
  }
  if (rows==columns)
    rows=rows/2;
  else
    columns=rows;
}

nautyVerify(graphName,0.0,0,0,0,0,0,0);
deallocate();
}

void hypercube(k)
int k;
{
int i,j;
double autos;
char string1[20],string2[20];

sprintf(graphName,"hypercube(%d)",k);
numVertices=twoPower(k);
allocate();
for (i=0;i<numVertices;i++)
  for (j=0;j<k;j++)
  {
    sprintf(string1,"%d",i);
    sprintf(string2,"%d",i^(twoPower(j)));
    connectVertices(string1,string2);
  }
autos=twoPower(k);
for (i=0;i<k;i++)
  autos*=(double)(i+1);
nautyVerify(graphName,autos,0,1,k,k,k,twoPower(k-1));
deallocate();
}

void shuffleExchange(k)
int k;
{
int i,j;
char string1[20],string2[20];

sprintf(graphName,"shuffleExchange(%d)",k);
numVertices=twoPower(k);
allocate();
for (i=0;i<numVertices-1;i++)
{
  sprintf(string1,"%d",i);
  sprintf(string2,"%d",i^1);
  connectVertices(string1,string2);
  sprintf(string2,"%d",2*i%(numVertices-1));
  connectVertices(string1,string2);
}
sprintf(string1,"%d",numVertices-1);
sprintf(string2,"%d",(numVertices-1)^1);
connectVertices(string1,string2);
connectVertices(string1,string1);
if (k==0)
  nautyVerify(graphName,1.0,0,1,1,2,2,0);
else
  nautyVerify(graphName,2.0,0,numVertices/2,2*k-1,2,2,0);
deallocate();
}

void undirectedShuffleExchange(k)
int k;
{
int i,j;
char string1[20],string2[20];

sprintf(graphName,"undirectedShuffleExchange(%d)",k);
numVertices=twoPower(k);
allocate();
for (i=0;i<numVertices-1;i++)
{
  sprintf(string1,"%d",i);
  sprintf(string2,"%d",i^1);
  connectVertices(string1,string2);
  sprintf(string2,"%d",2*i%(numVertices-1));
  if (strcmp(string1,string2)!=0)
    connectVertices(string1,string2);
/* Make shuffle links bi-directional to help METIS */
  if (strcmp(string1,string2)!=0)
    connectVertices(string2,string1);
}
sprintf(string1,"%d",numVertices-1);
sprintf(string2,"%d",(numVertices-1)^1);
connectVertices(string1,string2);
if (k==0)
  nautyVerify(graphName,1.0,0,1,1,2,2,0);
else if (k==1)
  nautyVerify(graphName,2.0,0,1,1,1,1,1);
else if (k==2)
  nautyVerify(graphName,2.0,0,2,3,2,2,1);
else
  nautyVerify(graphName,4.0,0,0,2*k-1,3,3,0);
deallocate();
}

void generalShuffleExchange(k)
int k;
{
int i,j;
char string1[20],string2[20];

sprintf(graphName,"generalShuffleExchange(%d)",k);
numVertices=k;
allocate();
for (i=0;i<numVertices/2;i++)
{
  sprintf(string1,"%d",i);
  if (i>0)
  {
    sprintf(string2,"%d",i-1);
    connectVertices(string1,string2);
  }
  sprintf(string2,"%d",i+1);
  connectVertices(string1,string2);
  sprintf(string2,"%d",2*i);
  connectVertices(string1,string2);
}
for (i=numVertices/2;i<numVertices;i++)
{
  sprintf(string1,"%d",i);
  sprintf(string2,"%d",i-1);
  connectVertices(string1,string2);
  if (i<numVertices-1)
  {
    sprintf(string2,"%d",i+1);
    connectVertices(string1,string2);
  }
  sprintf(string2,"%d",2*i-(numVertices-1));
  connectVertices(string1,string2);
}
if (k<6)
  nautyVerify(graphName,2.0,0,numVertices/2,0,2,2,0);
else
  nautyVerify(graphName,2.0,0,numVertices/2,0,3,3,0);
deallocate();
}

void deBruijn(k)
int k;
{
int i,j;
char string1[20],string2[20];

sprintf(graphName,"deBruijn(%d)",k);
numVertices=twoPower(k);
allocate();
for (i=0;i<numVertices;i++)
{
  sprintf(string1,"%d",i);
  sprintf(string2,"%d",(2*i)%numVertices);
  connectVertices(string1,string2);
  sprintf(string2,"%d",(2*i+1)%numVertices);
  connectVertices(string1,string2);
}
if (k==0)
  nautyVerify(graphName,1.0,0,(numVertices+1)/2,0,1,1,0);
else
  nautyVerify(graphName,2.0,0,(numVertices+1)/2,k,2,2,0);
deallocate();
}

void pm2i(k)
int k;
{
int i,j;
double autos;
char string1[20],string2[20];

sprintf(graphName,"pm2i(%d)",k);
numVertices=twoPower(k);
allocate();
for (j=0;j<numVertices;j++)
  for (i=0;i<k;i++)
  {
    sprintf(string1,"%d",j);
    sprintf(string2,"%d",(j+twoPower(i))%numVertices);
    connectVertices(string1,string2);
    sprintf(string2,"%d",(j-twoPower(i)+numVertices)%numVertices);
    connectVertices(string1,string2);
  }
if (k==0)
  autos=1.0;
else if (k==1)
  autos=2.0;
else if (k==2)
  autos=24.0;
else
  autos=(double)(twoPower(k+1));
nautyVerify(graphName,autos,0,1,(k+1)/2,2*k-1,2*k-1,twoPower(k));
deallocate();
}

void twoDmesh(n)
int n;
{
int i,j;
char string1[20],string2[20];

sprintf(graphName,"twoDmesh(%d)",n);
numVertices=n*n;
allocate();
for (i=0;i<n;i++)
  for (j=0;j<n;j++)
  {
    sprintf(string1,"row%dcol%d",i,j);
    if (i-1>=0)
    {
      sprintf(string2,"row%dcol%d",i-1,j);
      connectVertices(string1,string2);
    }
    if (i+1<n)
    {
      sprintf(string2,"row%dcol%d",i+1,j);
      connectVertices(string1,string2);
    }
    if (j-1>=0)
    {
      sprintf(string2,"row%dcol%d",i,j-1);
      connectVertices(string1,string2);
    }
    if (j+1<n)
    {
      sprintf(string2,"row%dcol%d",i,j+1);
      connectVertices(string1,string2);
    }
  }
if (n==1)
  nautyVerify(graphName,1.0,0,1,0,0,0,0);
else if (n==2)
  nautyVerify(graphName,8.0,0,1,2,2,2,2);
else if (n%2==0)
  nautyVerify(graphName,8.0,0,((n+1)/2)*((n+1)/2+1)/2,
    2*(n-1),4,4,n);
else
  nautyVerify(graphName,8.0,0,((n+1)/2)*((n+1)/2+1)/2,
    2*(n-1),4,4,n+1);
deallocate();
}

void twoDmesh2(m,n)
int m,n;
{
int i,j,vec;
double autos;
int diameter,indegree,outdegree,bisectionWidth;
char string1[20],string2[20];

if (m>n)
  abort();
sprintf(graphName,"twoDmesh2(%d,%d)",m,n);
numVertices=m*n;
allocate();
for (i=0;i<m;i++)
  for (j=0;j<n;j++)
  {
    sprintf(string1,"row%dcol%d",i,j);
    if (i-1>=0)
    {
      sprintf(string2,"row%dcol%d",i-1,j);
      connectVertices(string1,string2);
    }
    if (i+1<m)
    {
      sprintf(string2,"row%dcol%d",i+1,j);
      connectVertices(string1,string2);
    }
    if (j-1>=0)
    {
      sprintf(string2,"row%dcol%d",i,j-1);
      connectVertices(string1,string2);
    }
    if (j+1<n)
    {
      sprintf(string2,"row%dcol%d",i,j+1);
      connectVertices(string1,string2);
    }
  }
if (m==n)
  if (m==1)
  {
    autos=1.0;
    vec=1;
    diameter=0;
    indegree=outdegree=0;
    bisectionWidth=0;
  }
  else
  {
    autos=8.0;
    vec=((n+1)/2)*((n+1)/2+1)/2;
    diameter=2*(m-1);
    if (m==2)
      indegree=outdegree=2;
    else
      indegree=outdegree=4;
    if (m%2==0)
      bisectionWidth=m;
    else
      bisectionWidth=m+1;
  }
else
{
  diameter=m+n-2;
  if (m==1)
  {
    autos=2.0;
    if (n==2)
      indegree=outdegree=1;
    else
      indegree=outdegree=2;
    bisectionWidth=1;
  }
  else
  {
    autos=4.0;
    if (m==2)
      indegree=outdegree=3;
    else
      indegree=outdegree=4;
    if (n%2==0)
      bisectionWidth=m;
    else
      bisectionWidth=m+1;
  }
  vec=((m+1)/2)*((n+1)/2);
}
nautyVerify(graphName,autos,0,vec,diameter,indegree,outdegree,
  bisectionWidth);
deallocate();
}

void twoDtorus(n)
int n;
{
int i,j;
double autos;
int diameter,indegree,outdegree,bisectionWidth;
char string1[20],string2[20];

sprintf(graphName,"twoDtorus(%d)",n);
numVertices=n*n;
allocate();
for (i=0;i<n;i++)
  for (j=0;j<n;j++)
  {
    sprintf(string1,"row%dcol%d",i,j);
    sprintf(string2,"row%dcol%d",(i-1+n)%n,j);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",(i+1)%n,j);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",i,(j-1+n)%n);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",i,(j+1)%n);
    connectVertices(string1,string2);
  }
if (n==1)
{
  diameter=indegree=outdegree=bisectionWidth=0;
  autos=1.0;
}
else if (n==2)
{
  diameter=2;
  indegree=outdegree=bisectionWidth=2;
  autos=8.0;
}
else
{
  if (n%2==0)
    diameter=n;
  else
    diameter=n-1;
  indegree=outdegree=4;
  if (n%2==0)
    bisectionWidth=2*n;
  else
    bisectionWidth=2*(n+1);
  if (n==4)
    autos=384.0;
  else
    autos=8*n*n;
}
nautyVerify(graphName,autos,0,1,diameter,indegree,outdegree,bisectionWidth);
deallocate();
}

void twoDtorus2(m,n)
int m,n;
{
int i,j;
double autos;
int diameter,indegree,outdegree,bisectionWidth;
char string1[20],string2[20];

if (m>n)
  abort();
sprintf(graphName,"twoDtorus2(%d,%d)",m,n);
numVertices=m*n;
allocate();
for (i=0;i<m;i++)
  for (j=0;j<n;j++)
  {
    sprintf(string1,"row%dcol%d",i,j);
    sprintf(string2,"row%dcol%d",(i-1+m)%m,j);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",(i+1)%m,j);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",i,(j-1+n)%n);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",i,(j+1)%n);
    connectVertices(string1,string2);
  }
if (m==n)
{
  if (n==1)
  {
    diameter=indegree=outdegree=bisectionWidth=0;
    autos=1.0;
  }
  else if (n==2)
  {
    diameter=2;
    indegree=outdegree=bisectionWidth=2;
    autos=8.0;
  }
  else
  {
    if (n%2==0)
      diameter=n;
    else
      diameter=n-1;
    indegree=outdegree=4;
    if (n%2==0)
      bisectionWidth=2*n;
    else
      bisectionWidth=2*(n+1);

    if (n==4)
      autos=384.0;
    else
      autos=8*n*n;
  }
}
else if (m==1)
  if (n==2)
  {
    diameter=1;
    indegree=outdegree=2;
    bisectionWidth=2;
    autos=2.0;
  }
  else
  {
    autos=2.0*n;
    diameter=n/2;
    indegree=outdegree=3;
    bisectionWidth=2;
  }
else
{
  diameter=m/2+n/2;
  if (m==2)
    indegree=outdegree=3;
  else
    indegree=outdegree=4;
  if (n%2==0)
    bisectionWidth=2*m;
  else
    bisectionWidth=2*(m+1);
  if (m==2)
  {
    if (n==3)
      autos=12.0;
    else if (n==4)
      autos=48.0;
    else
      autos=4*n;
  }
  else
    autos=4*m*n;
}
nautyVerify(graphName,autos,0,1,diameter,indegree,outdegree,bisectionWidth);
deallocate();
}

void linear(n)
int n;
{
int i,j;
char string1[20],string2[20];

sprintf(graphName,"linear(%d)",n);
numVertices=n;
allocate();
for (i=0;i<n-1;i++)
{
  sprintf(string1,"%d",i);
  sprintf(string2,"%d",i+1);
  connectVertices(string1,string2);
  connectVertices(string2,string1);
}
if (n==1)
  nautyVerify(graphName,1.0,0,1,0,0,0,0);
else if (n==2)
  nautyVerify(graphName,2.0,0,1,1,1,1,1);
else
  nautyVerify(graphName,2.0,0,(n+1)/2,n-1,2,2,1);
deallocate();
}

void ring(n)
int n;
{
int i,j;
char string1[20],string2[20];

sprintf(graphName,"ring(%d)",n);
numVertices=n;
allocate();
for (i=0;i<n;i++)
{
  sprintf(string1,"%d",i);
  sprintf(string2,"%d",(i+1)%n);
  connectVertices(string1,string2);
  connectVertices(string2,string1);
}
if (n==1)
  nautyVerify(graphName,1.0,0,1,0,0,0,0);
else if (n==2)
  nautyVerify(graphName,2.0,0,1,1,1,1,1);
else
  nautyVerify(graphName,2.0*n,0,1,n/2,2,2,2);
deallocate();
}

void wrappedButterfly(numColumns)
int numColumns;
{
int numRows,i,j;
double autos;
int diameter,indegree,outdegree;
char string1[20],string2[20];

sprintf(graphName,"wrappedButterfly(%d)",numColumns);
numRows=twoPower(numColumns);
numVertices=numRows*numColumns;
allocate();
for (i=0;i<numRows;i++)
  for (j=0;j<numColumns;j++)
  {
    sprintf(string1,"row%dcol%d",i,j);
    sprintf(string2,"row%dcol%d",i,(j+1)%numColumns);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",i,(j-1+numColumns)%numColumns);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",i^twoPower(j),(j+1)%numColumns);
    connectVertices(string1,string2);
    connectVertices(string2,string1);
  }
diameter=numColumns+numColumns/2;
if (numColumns==1)
{
  autos=2.0;
  indegree=outdegree=2;
}
else if (numColumns==2)
{
  autos=48.0;
  indegree=outdegree=3;
}
else
{
  autos=numColumns*(twoPower(numColumns+1));
  indegree=outdegree=4;
}
nautyVerify(graphName,autos,0,1,diameter,indegree,outdegree,
  twoPower(numColumns));
deallocate();
}

void cubeConnection(numColumns)
int numColumns;
{
int numRows,i,j;
char string1[20],string2[20];

sprintf(graphName,"cubeConnection(%d)",numColumns);
numRows=twoPower(numColumns);
numVertices=numRows*numColumns;
allocate();
for (i=0;i<numRows;i++)
  for (j=0;j<numColumns;j++)
  {
    sprintf(string1,"row%dcol%d",i,j);
    sprintf(string2,"row%dcol%d",i,(j+1)%numColumns);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",i,(j-1+numColumns)%numColumns);
    connectVertices(string1,string2);
    sprintf(string2,"row%dcol%d",i^twoPower(j),j);
    connectVertices(string1,string2);
    connectVertices(string2,string1);
  }
if (numColumns==1)
  nautyVerify(graphName,2.0,0,1,0,0,0,0);
else
  nautyVerify(graphName,(double)numColumns*(twoPower(numColumns+1)),0,1,
    2*numColumns,3,3,twoPower(numColumns-1));
deallocate();
}

void butterfly(numColumns)
int numColumns;
{
int numRows,i,j;
char string1[20],string2[20];

sprintf(graphName,"butterfly(%d)",numColumns);
numRows=twoPower(numColumns-1);
numVertices=numRows*numColumns;
allocate();
for (i=0;i<numRows;i++)
  for (j=0;j<numColumns-1;j++)
  {
    sprintf(string1,"row%dcol%d",i,j);
    sprintf(string2,"row%dcol%d",i,j+1);
    connectVertices(string1,string2);
    connectVertices(string2,string1);
    sprintf(string2,"row%dcol%d",i^twoPower(j),j+1);
    connectVertices(string1,string2);
    connectVertices(string2,string1);
  }
if (numColumns==1)
  nautyVerify(graphName,1.0,0,1,0,0,0,0);
else if (numColumns==2)
  nautyVerify(graphName,8.0,0,1,2,2,2,2);
else
  nautyVerify(graphName,twoPowerBig((twoPower(numColumns))-1),
    0,(numColumns+1)/2,2*(numColumns-1),4,4,twoPower(numColumns-1));
deallocate();
}

/* The following code was written by Kenneth Damrau, 2/98 */
/* Evaluates k-partite graphs analytically and by search */

/*
 * The first argument is the number of sets of vertices (of size x_i)
 * with an edge joining two vertices if and only if they lie in
 * different sets.  The rest of the arguments are the sizes x_i.
 *
 * called like k_partite(2, 4, 5)
 * or     like k_partite(5, 2, 3, 5, 3, 2);
 */

int compar (const void *p, const void *q);
int factorial (int n);
int tothe (int x, int y);       /* returns x to the y */

void
k_partite (int k,...)
{
  int BeginningOfSet, EndOfSet;
  int i, j, SetNumber;
  int *s, *p, *x;
  int product = 1;
  int Y;                        /* number of different sizes */
  va_list q;
  char *ptr;
  int i1,j1;
  char string1[25],string2[25];

  sprintf(graphName,"k_partite(%d,",k);
  ptr=graphName+strlen(graphName);
  numVertices = 0;
  s = calloc (k, sizeof (int));
  p = calloc (k, sizeof (int));
  x = malloc (k * sizeof (int));
  if (s == NULL || p == NULL || x == NULL)
    {
      fprintf (stderr, "malloc failed file %s line %d\n", __FILE__, __LINE__);
      exit (EXIT_FAILURE);
    }
  va_start (q, k);
  for (i = 0; i < k; i++)
    {
      x[i] = va_arg (q, int);
      numVertices += x[i];
      if (i<k-1)
      {
        sprintf(ptr,"%d,",x[i]);
        ptr+=strlen(ptr);
      }
      else
        sprintf(ptr,"%d)",x[i]);
    }
  va_end (q);
  qsort (x, k, sizeof (int), compar);
  BeginningOfSet = 0;
  EndOfSet = x[0] - 1;
  allocate ();
  for (SetNumber = 0; SetNumber < k; SetNumber++)
    {
      for (i = BeginningOfSet; i <= EndOfSet; i++)
        {
          for (j = EndOfSet + 1; j < numVertices; j++)
            {
/*
              matrix[i][j] = 1;
              matrix[j][i] = 1;
*/
            }
        }
      BeginningOfSet = EndOfSet + 1;
      EndOfSet += x[SetNumber + 1];
    }
/* Code for labeled vertices */
  for (i=0;i<k;i++)
    for (j=i+1;j<k;j++)
      for (i1=0;i1<x[i];i1++)
        for (j1=0;j1<x[j];j1++)
        {
          sprintf(string1,"%d.%d",i,i1);
          sprintf(string2,"%d.%d",j,j1);
          connectVertices(string1,string2);
          connectVertices(string2,string1);
        }

  for (i = 0; i < k; i++)
    {
      if (i == 0)
        {
          Y = 1;
          s[0] = x[0];
        }
      else if (x[i] != x[i - 1])
        {
          s[Y] = x[i];
          ++Y;
        }
      ++p[Y - 1];
    }
  for (i = 0; i < Y; i++)
    product *= factorial (p[i]) * tothe (factorial (s[i]), p[i]);
/*
  printf ("automorphisms = %d classes = %d\n", product, Y);
*/
  nautyVerify(graphName,(double) product,0,Y,0,0,0,0);
  deallocate();
  free (s);
  free (p);
  free (x);
  fflush(stdout);
}

int
factorial (int n)
{
  if (n == 0 || n == 1)
    return 1;
  if (n == 2)
    return 2;
  return (n * factorial (n - 1));
}


int
compar (const void *p, const void *q)
{
  return *((int *) p) - *((int *) q);
}

int
tothe (int x, int y)
{
  int i, product = 1;

  for (i = 0; i < y; ++i)
    product *= x;
  return product;
}

/* End of k-partite code */

/* Code for ranking and unranking permutations */

long long columnValue[50];

long long cardPermset(n,k)
int n,k;
{
int i;
long long card;

columnValue[k-1]=1;
for (i=k-2;i>=0;i--)
  columnValue[i]=columnValue[i+1]*(n-i-1);
card=n*columnValue[0];
/*
printf("column values:\n");
for (i=0;i<k;i++)
  printf("%f\n",(double) columnValue[i]);
printf("cardinality of permutation set %f\n",(double) card);
*/
return card;
}

long long perm2rank(perm,n,k)
int *perm,n,k;
{
int i,j;
long long rank;
int available[50],r[50];

for (i=0;i<n;i++)
  available[i]=i;
for (i=0;i<k;i++)
{
  for (j=0;perm[i]!=available[j];j++)
    ;
  r[i]=j;
  for (;j<n-1-i;j++)
    available[j]=available[j+1];
}
rank=0;
for (i=0;i<k;i++)
  rank+=r[i]*columnValue[i];
return rank;
}

void rank2perm(rank,perm,n,k)
long long rank;
int *perm,n,k;
{
int i,j;
int r[50],available[50];

for (i=0;i<k;i++)
{
  r[i]=rank/columnValue[i];
  rank=rank%columnValue[i];
}

for (i=0;i<n;i++)
  available[i]=i;
for (i=0;i<k;i++)
{
  perm[i]=available[r[i]];
  for (j=r[i];j<n-1-i;j++)
    available[j]=available[j+1];
}
}

void processGroupGraph(int numGenerators,int permSize,int **generator)
{
long long hashRank[8000],rank,newRank;
int i,j,k;
long long queue[8000];
int tail=(-1),head=(-1);
int workPerm[50],newPerm[50];
int tailVertex,headVertex;
char tailString[500],headString[500],*ptr;
float startCPU;

startCPU=CPUseconds();
numVertices=0;
for (i=0;i<8000;i++)
  hashRank[i]=(-1);
cardPermset(permSize,permSize);
for (i=0;i<numGenerators;i++)
{
  rank=perm2rank(generator[i],permSize,permSize);
  j=rank%8000;
  while (hashRank[j]!=(-1) && hashRank[j]!=rank)
    j=(j+1)%8000;
  if (hashRank[j]==rank)
    continue;
  hashRank[j]=rank;
  numVertices++;
  tail++;
  queue[tail]=rank;
}
while (tail!=head)
{
  head++;
  rank2perm(queue[head],workPerm,permSize,permSize);
  for (i=0;i<numGenerators;i++)
  {
    for (j=0;j<permSize;j++)
      newPerm[j]=generator[i][workPerm[j]];
    newRank=perm2rank(newPerm,permSize,permSize);
    j=newRank%8000;
    while (hashRank[j]!=(-1) && hashRank[j]!=newRank)
      j=(j+1)%8000;
    if (hashRank[j]==newRank)
      continue;
    hashRank[j]=newRank;
    numVertices++;
    tail++;
    queue[tail]=newRank;
  }
}
permutationCPU+=CPUseconds()-startCPU;
allocate();
startCPU=CPUseconds();
for (i=0;i<8000;i++)
  if (hashRank[i]!=(-1))
  {
    rank2perm(hashRank[i],workPerm,permSize,permSize);
    sprintf(tailString,"perm%d",workPerm[0]);
    ptr=tailString+strlen(tailString);
    for (j=1;j<permSize;j++)
    {
      sprintf(ptr,".%d",workPerm[j]);
      ptr+=strlen(ptr);
    }
    for (j=0;j<numGenerators;j++)
    {
      for (k=0;k<permSize;k++)
        newPerm[k]=generator[j][workPerm[k]];
      sprintf(headString,"perm%d",newPerm[0]);
      ptr=headString+strlen(headString);
      for (k=1;k<permSize;k++)
      {
        sprintf(ptr,".%d",newPerm[k]);
        ptr+=strlen(ptr);
      }
      permutationCPU+=CPUseconds()-startCPU;
      connectVertices(tailString,headString);
      startCPU=CPUseconds();
    }
  }
}

void starGraph(int k)
{
int **generator,i,j;
double autos;

sprintf(graphName,"starGraph(%d)",k);
generator=(int**) malloc((k-1)*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<k-1;i++)
  generator[i]=(int*) malloc(k*sizeof(int));
for (i=0;i<k-1;i++)
{
  for (j=0;j<k;j++)
    generator[i][j]=j;
  generator[i][0]=i+1;
  generator[i][i+1]=0;
}
processGroupGraph(k-1,k,generator);
for (i=0;i<k-1;i++)
  free(generator[i]);
free(generator);
autos=factorial(k)*factorial(k-1);
nautyVerify(graphName,autos,0,1,k,k-1,k-1,factorial(k)/3);
deallocate();
}

void bubbleSortGraph(int k)
{
int **generator,i,j;
double autos;

sprintf(graphName,"bubbleSortGraph(%d)",k);
generator=(int**) malloc((k-1)*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<k-1;i++)
  generator[i]=(int*) malloc(k*sizeof(int));
for (i=0;i<k-1;i++)
{
  for (j=0;j<k;j++)
    generator[i][j]=j;
  generator[i][i]=i+1;
  generator[i][i+1]=i;
}
processGroupGraph(k-1,k,generator);
for (i=0;i<k-1;i++)
  free(generator[i]);
free(generator);
if (k==2)
  autos=2.0;
else
  autos=2.0*factorial(k);
nautyVerify(graphName,autos,0,1,k*(k-1)/2,k-1,k-1,0);
deallocate();
}

void modifiedBubbleSortGraph(int k)
{
int **generator,i,j;
double autos;

sprintf(graphName,"modifiedBubbleSortGraph(%d)",k);
generator=(int**) malloc(k*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<k;i++)
  generator[i]=(int*) malloc(k*sizeof(int));
for (i=0;i<k-1;i++)
{
  for (j=0;j<k;j++)
    generator[i][j]=j;
  generator[i][i]=i+1;
  generator[i][i+1]=i;
}
for (j=0;j<k;j++)
  generator[k-1][j]=j;
generator[k-1][0]=k-1;
generator[k-1][k-1]=0;
processGroupGraph(k,k,generator);
for (i=0;i<k;i++)
  free(generator[i]);
free(generator);
if (k==2)
  autos=2.0;
else if (k==3)
  autos=72.0;
else if (k==4)
  autos=768.0;
else
  autos=2.0*k*factorial(k);
nautyVerify(graphName,autos,0,1,0,k,k,2*factorial(k-1));
deallocate();
}

void prefixReversal(int k)
{
int **generator,i,j;
double autos;

sprintf(graphName,"prefixReversal(%d)",k);
generator=(int**) malloc((k-1)*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<k-1;i++)
  generator[i]=(int*) malloc(k*sizeof(int));
for (i=0;i<k-1;i++)
{
  for (j=0;j<k;j++)
    generator[i][j]=j;
  for (j=0;j<=i+1;j++)
    generator[i][j]=i+1-j;
}
processGroupGraph(k-1,k,generator);
for (i=0;i<k-1;i++)
  free(generator[i]);
free(generator);
if (k==3 || k==4)
  autos=2.0*factorial(k);
else
  autos=factorial(k);
nautyVerify(graphName,autos,0,1,0,k-1,k-1,factorial(k)/3);
deallocate();
}

void completeTransposition(int k)
{
int **generator,i,j,l,m,numGenerators;
double autos;

sprintf(graphName,"completeTransposition(%d)",k);
numGenerators=k*(k-1)/2;
generator=(int**) malloc(numGenerators*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<numGenerators;i++)
  generator[i]=(int*) malloc(k*sizeof(int));
l=0;
for (i=0;i<k-1;i++)
  for (j=i+1;j<k;j++)
  {
    for (m=0;m<k;m++)
      generator[l][m]=m;
    generator[l][i]=j;
    generator[l][j]=i;
    l++;
  }
processGroupGraph(numGenerators,k,generator);
for (i=0;i<numGenerators;i++)
  free(generator[i]);
free(generator);
if (k==2)
  autos=2.0;
else
  autos=2.0*factorial(k)*factorial(k);
nautyVerify(graphName,autos,0,1,k-1,k*(k-1)/2,k*(k-1)/2,0);
deallocate();
}

void hypercubeGroup(k)
int k;
{
int i,j;
int **generator;
double autos;

sprintf(graphName,"hypercubeGroup(%d)",k);
generator=(int**) malloc(k*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<k;i++)
  generator[i]=(int*) malloc(2*k*sizeof(int));
for (i=0;i<k;i++)
{
  for (j=0;j<2*k;j++)
    generator[i][j]=j;
  generator[i][2*i]=2*i+1;
  generator[i][2*i+1]=2*i;
}
processGroupGraph(k,2*k,generator);
for (i=0;i<k;i++)
  free(generator[i]);
free(generator);
autos=twoPowerBig(k);
for (i=0;i<k;i++)
  autos*=(double)(i+1);
nautyVerify(graphName,autos,0,1,0,0,0,0);
deallocate();
}

void ringGroup(n)
int n;
{
int i,j;
int **generator;

sprintf(graphName,"ringGroup(%d)",n);
generator=(int**) malloc(2*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<2;i++)
  generator[i]=(int*) malloc(n*sizeof(int));
for (i=0;i<n;i++)
  generator[0][i]=(i+1)%n;
for (i=0;i<n;i++)
  generator[1][i]=(i-1+n)%n;
processGroupGraph(2,n,generator);
for (i=0;i<2;i++)
  free(generator[i]);
free(generator);
if (n==1)
  nautyVerify(graphName,1.0,0,1,0,0,0,0);
else if (n==2)
  nautyVerify(graphName,2.0,0,1,0,0,0,0);
else
  nautyVerify(graphName,2.0*n,0,1,0,0,0,0);
deallocate();
}

void twoDtorusGroup(m,n)
int m,n;
{
int i,j;
int **generator;
double autos;

if (m>n)
  abort();
sprintf(graphName,"twoDtorusGroup(%d,%d)",m,n);
generator=(int**) malloc(4*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<4;i++)
  generator[i]=(int*) malloc((m+n)*sizeof(int*));
for (i=0;i<4;i++)
  for (j=0;j<m+n;j++)
    generator[i][j]=j;
for (i=0;i<n;i++)
  generator[0][i]=(i+1)%n;
for (i=0;i<n;i++)
  generator[1][i]=(i-1+n)%n;
for (i=0;i<m;i++)
  generator[2][n+i]=n+(i+1)%m;
for (i=0;i<m;i++)
  generator[3][n+i]=n+(i-1+m)%m;
processGroupGraph(4,m+n,generator);
for (i=0;i<4;i++)
  free(generator[i]);
free(generator);
if (m==n)
{
  if (n==1)
    autos=1.0;
  else if (n==2)
    autos=8.0;
  else if (n==4)
    autos=384.0;
  else
    autos=8*n*n;
}
else if (m==1)
  if (n==2)
    autos=2;
  else
    autos=2*n;
else if (m==2)
{
  if (n==3)
    autos=12.0;
  else if (n==4)
    autos=48.0;
  else
    autos=4*n;
}
else
  autos=4*m*n;
nautyVerify(graphName,autos,0,1,0,0,0,0);
deallocate();
}

void threeDtorusGroup(m,n,p)
int m,n,p;
{
int i,j;
int **generator;
double autos;

if (m>n || n>p)
  abort();
sprintf(graphName,"threeDtorusGroup(%d,%d,%d)",m,n,p);
generator=(int**) malloc(6*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<6;i++)
  generator[i]=(int*) malloc((m+n+p)*sizeof(int*));
for (i=0;i<6;i++)
  for (j=0;j<m+n+p;j++)
    generator[i][j]=j;
for (i=0;i<n;i++)
  generator[0][i]=(i+1)%n;
for (i=0;i<n;i++)
  generator[1][i]=(i-1+n)%n;
for (i=0;i<m;i++)
  generator[2][n+i]=n+(i+1)%m;
for (i=0;i<m;i++)
  generator[3][n+i]=n+(i-1+m)%m;
for (i=0;i<p;i++)
  generator[4][n+m+i]=n+m+(i+1)%p;
for (i=0;i<p;i++)
  generator[5][n+m+i]=n+m+(i-1+p)%p;
processGroupGraph(6,m+n+p,generator);
for (i=0;i<6;i++)
  free(generator[i]);
free(generator);
nautyVerify(graphName,0.0,0,1,0,0,0,0);
deallocate();
}

void cubeConnectionGroup(numColumns)
int numColumns;
{
int numRows,i,j;
int *generator[3];

sprintf(graphName,"cubeConnectionGroup(%d)",numColumns);
for (i=0;i<3;i++)
  generator[i]=(int*) malloc(2*numColumns*sizeof(int));
for (i=0;i<numColumns;i++)
  generator[0][i]=(i+1)%numColumns;
for (i=0;i<numColumns;i++)
  generator[0][numColumns+i]=numColumns+(i+1)%numColumns;
for (i=0;i<numColumns;i++)
  generator[1][i]=(i-1+numColumns)%numColumns;
for (i=0;i<numColumns;i++)
  generator[1][numColumns+i]=numColumns+(i-1+numColumns)%numColumns;
for (i=0;i<2*numColumns;i++)
  generator[2][i]=i;
generator[2][0]=numColumns;
generator[2][numColumns]=0;
processGroupGraph(3,2*numColumns,generator);
for (i=0;i<3;i++)
  free(generator[i]);
if (numColumns==1)
  nautyVerify(graphName,2.0,0,1,0,0,0,0);
else
  nautyVerify(graphName,(double)numColumns*(twoPower(numColumns+1)),0,1,
    0,0,0,0);
deallocate();
}

void wrappedButterflyGroup(numColumns)
int numColumns;
{
int numRows,i,j;
double autos;
int *generator[6],temp;

sprintf(graphName,"wrappedButterflyGroup(%d)",numColumns);
for (i=0;i<6;i++)
  generator[i]=(int*) malloc(2*numColumns*sizeof(int));
for (i=0;i<numColumns;i++)
  generator[0][i]=(i+1)%numColumns;
for (i=0;i<numColumns;i++)
  generator[0][numColumns+i]=numColumns+(i+1)%numColumns;
for (i=0;i<numColumns;i++)
  generator[1][i]=(i-1+numColumns)%numColumns;
for (i=0;i<numColumns;i++)
  generator[1][numColumns+i]=numColumns+(i-1+numColumns)%numColumns;
for (i=0;i<2*numColumns;i++)
  generator[2][i]=generator[0][i];
temp=generator[2][0];
generator[2][0]=generator[2][numColumns];
generator[2][numColumns]=temp;
for (i=0;i<2*numColumns;i++)
  generator[3][i]=generator[0][i];
temp=generator[3][1];
generator[3][1]=generator[3][numColumns+1];
generator[3][numColumns+1]=temp;
for (i=0;i<2*numColumns;i++)
  generator[4][i]=generator[1][i];
temp=generator[4][0];
generator[4][0]=generator[4][numColumns];
generator[4][numColumns]=temp;
for (i=0;i<2*numColumns;i++)
  generator[5][i]=generator[1][i];
temp=generator[5][1];
generator[5][1]=generator[5][numColumns+1];
generator[5][numColumns+1]=temp;
processGroupGraph(6,2*numColumns,generator);
for (i=0;i<6;i++)
  free(generator[i]);
if (numColumns==1)
  autos=2.0;
else if (numColumns==2)
  autos=48.0;
else
  autos=numColumns*(twoPower(numColumns+1));
nautyVerify(graphName,autos,0,1,0,0,0,0);
deallocate();
}

void pm2iGroup(k)
int k;
{
int i,j;
double autos;
int **generator;

sprintf(graphName,"pm2iGroup(%d)",k);
generator=(int**) malloc((2*k-1)*sizeof(int*));
if (!generator)
  memoryDie();
for (i=0;i<2*k-1;i++)
  generator[i]=(int*) malloc(twoPower(k)*sizeof(int));
for (i=0;i<k;i++)
  for (j=0;j<twoPower(k);j++)
    generator[i][j]=(twoPower(i)+j)%twoPower(k);
for (i=0;i<k-1;i++)
  for (j=0;j<twoPower(k);j++)
    generator[k+i][j]=(twoPower(k)-twoPower(i)+j)%twoPower(k);
processGroupGraph(2*k-1,twoPower(k),generator);
for (i=0;i<2*k-1;i++)
  free(generator[i]);
free(generator);
if (k==0)
  autos=1.0;
else if (k==1)
  autos=2.0;
else if (k==2)
  autos=24.0;
else
  autos=(double)(twoPower(k+1));
nautyVerify(graphName,autos,0,1,0,0,0,0);
deallocate();
}

main()
{
int i,j,k,l;
float startCPU,endCPU;

startCPU=CPUseconds();

if ((scriptDescript=fopen("format.script","w"))==NULL)
  printf("format.script file open failed\n");

for (i=1;i<5;i++)
  CBTmesh(i);

/*
for (i=1;i<=127;i++)
  tree(i);
for (i=1;i<=6;i++)
  for (j=0;j<=6;j++)
    linearCBT(i,j);
for (i=0;i<=4;i++)
  for (j=i;j<=4;j++)
    meshCBT(i,j);
for (i=1;i<=12;i++)
  hypercube(i);
for (i=1;i<=10;i++)
  shuffleExchange(i);
for (i=1;i<=10;i++)
  undirectedShuffleExchange(i);
for (i=2;i<=64;i+=2)
  generalShuffleExchange(i);
for (i=0;i<=6;i++)
  deBruijn(i);
for (i=1;i<9;i++)
  pm2i(i);
for (i=2;i<=10;i++)
  twoDmesh(i);
for (i=1;i<=10;i++)
  for (j=i;j<=10;j++)
    if (i!=1 || j!=1)
      twoDmesh2(i,j);
for (i=2;i<=15;i++)
  twoDtorus(i);
for (i=1;i<=20;i++)
  for (j=i;j<=20;j++)
    if (i!=1 || j!=1)
      twoDtorus2(i,j);
for (i=2;i<=10;i++)
  linear(i);
for (i=2;i<=10;i++)
  ring(i);
for (i=1;i<=9;i++)
  wrappedButterfly(i);
for (i=1;i<=9;i++)
  cubeConnection(i);
for (i=1;i<=9;i++)
  butterfly(i);
for (i=1;i<=4;i++)
  for (j=i;j<=4;j++)
    k_partite(2,i,j);
for (i=1;i<=4;i++)
  for (j=i;j<=4;j++)
    for (k=j;k<=4;k++)
      k_partite(3,i,j,k);
for (i=1;i<=4;i++)
  for (j=i;j<=4;j++)
    for (k=j;k<=4;k++)
      for (l=k;l<=4;l++)
        k_partite(4,i,j,k,l);
for (i=2;i<=7;i++)
  starGraph(i);
for (i=2;i<=7;i++)
  bubbleSortGraph(i);
for (i=2;i<=7;i++)
  modifiedBubbleSortGraph(i);
for (i=2;i<=7;i++)
  prefixReversal(i);
for (i=2;i<=7;i++)
  completeTransposition(i);
for (i=1;i<=8;i++)
  hypercubeGroup(i);
for (i=2;i<=10;i++)
  ringGroup(i);
for (i=1;i<=10;i++)
  for (j=i;j<=10;j++)
    twoDtorusGroup(i,j);
for (i=1;i<=6;i++)
  for (j=i;j<=6;j++)
    for (k=j;k<=6;k++)
      threeDtorusGroup(i,j,k);
for (i=1;i<=7;i++)
  cubeConnectionGroup(i);
for (i=2;i<=4;i++)
  wrappedButterflyGroup(i);
for (i=1;i<=4;i++)
  pm2iGroup(i);
*/
endCPU=CPUseconds();
fclose(scriptDescript);
printf("CPU time: %f\n",endCPU-startCPU);
printf("nauty CPU: %f\n",nautyCPU);
printf("METIS CPU: %f\n",metisCPU);
printf("memory management CPU: %f\n",memoryCPU);
printf("hash CPU: %f\n",hashCPU);
printf("Lukes CPU: %f\n",LukesCPU);
printf("permutation CPU: %f\n",permutationCPU);
printf("qsort CPU: %f\n",sortCPU);
printf("unaccounted CPU: %f\n",endCPU-startCPU-nautyCPU-metisCPU
  -memoryCPU-hashCPU-LukesCPU-permutationCPU-sortCPU);
}