// Binary array multiplication using Kronrod's (4 russians)
// algorithm based on matrix rows partitioned into unsigned chars
// BPW, 6/15/07

// This code runs significantly faster if compiled with -O1 (at least).

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include <sys/resource.h>

#define SEED (666)

// Converts column bit subscript to subscript for unsigned char
#define SUB2CHAR(sub) ((sub)/8)

// Determines bit position within unsigned char for column bit subscript
#define SUB2BIT(sub) ((sub)%8)

// Gets 0/1 value for a column bit subscript
#define GETBIT(arr,row,col) ((arr[row][SUB2CHAR(col)]>>SUB2BIT(col))&1)

// Sets bit to 1
#define SETBIT(arr,row,col) \
  (arr[row][SUB2CHAR(col)]=arr[row][SUB2CHAR(col)]|(1<<SUB2BIT(col)))

// Clears bit to 0
#define CLEARBIT(arr,row,col) \
  (arr[row][SUB2CHAR(col)]=arr[row][SUB2CHAR(col)]&(255-(1<<SUB2BIT(col))))

float CPUtime()
{
  struct rusage rusage;
  getrusage(RUSAGE_SELF,&rusage);
  return rusage.ru_utime.tv_sec+rusage.ru_utime.tv_usec/1000000.0
    + rusage.ru_stime.tv_sec+rusage.ru_stime.tv_usec/1000000.0;
}

void charBinMM(int n,unsigned char a[256][256],unsigned char b[256][256],
  unsigned char c[256][256])
// Binary mat mult for n by n matrix a, n by n matrix b
// to produce n by n matrix c with each bit stored alone
// as an unsigned char.
{
int i,j,k;

n=256;

for (i=0;i<n;i++)
  for (j=0;j<n;j++)
  {
    c[i][j]=0;
    for (k=0;k<n;k++)
      if (a[i][k] && b[k][j])
      {
        c[i][j]=1;
        break;
      }
  }
}

void packedBinMM(int n,unsigned char A[256][32],unsigned char B[256][32],
  unsigned char C[256][32])
// Binary mat mult for n by n matrix A, n by n matrix B
// to produce n by n matrix C with 8 bits stored in each
// unsigned char.  (So the number of bytes in a row is n/8)
{
int i,j,k,p,t=8;

n=256;

for (i=0; i<n; i++)
{
  for (j=0; j<n/t; j++)       // Clear output row i
    C[i][j]=0;
  for (j=0; j<n/t; j++)
    for (k=0; k<t; k++)
      if ((A[i][j]>>k) & 1)   // Examine row i, bit # t*j+k
        for (p=0; p<n/t; p++) // Set corresponding bits
          C[i][p]=C[i][p] | B[t*j+k][p];
}
}

void KronrodBinMM(int n,unsigned char A[256][32],unsigned char B[256][32],
  unsigned char C[256][32])
// Binary mat mult for n by n matrix A, n by n matrix B
// to produce n by n matrix C with 8 bits stored in each
// unsigned char.  (So the number of bytes in each A row is n/8 and
// in each B row is p/8)

// Uses a byte-based version of Kronrod's algorithm.
{
int i,j,k,q,t=8;
// Kronrod subset rows
unsigned char subset[256][32];  // Work for combinations of or'd rows
int rangeStart;

n=256;

// Initialize output matrix
for (i=0; i<n; i++)
  for (j=0; j<n/t; j++)
    C[i][j]=0;

// Initialize empty subset - same for each group
for (i=0; i<n/t; i++)
  subset[0][i]=0;
// There are n/t groups of t rows in B
for (i=0; i<n/t; i++)
{
  rangeStart=1;
  // Include each row to achieve combinations OR'ed together
  for (j=0; j<t; j++)
  {
    for (k=rangeStart; k<2*rangeStart; k++)
      for (q=0; q<n/t; q++)
        // OR the subset without row j of this group with row j
        // of this group.  First row of group is B[t*i].
        subset[k][q]=subset[k-rangeStart][q] | B[t*i+j][q];
    rangeStart*=2;
  }

  // Update result rows based on this group of t B rows
  for (j=0; j<n; j++)
    for (q=0; q<n/t; q++)
      C[j][q]=C[j][q] | subset[A[j][i]][q];
}
}

main()
{
// Arrays for computing "unpacked" with one bit per unsigned char
unsigned char aSlow[256][256],bSlow[256][256],cSlow[256][256];
// Arrays for computing "packed" with 8 bits per unsigned char
unsigned char a[256][32],b[256][32],c[256][32];
int i,j,k,p;

float startCPU,stopCPU;

// Load conventional matrices with random bits
srandom(SEED);
for (i=0;i<256;i++)
  for (j=0;j<256;j++)
  {
    aSlow[i][j]=(0==abs(random())%50);
    bSlow[i][j]=(0==abs(random())%50);
  }

// Compute binary product
startCPU=CPUtime();
charBinMM(256,aSlow,bSlow,cSlow);
stopCPU=CPUtime();
printf("slow (unpacked) CPU %.5f\n",stopCPU-startCPU);

// Load packed matrices
srandom(SEED);
for (i=0;i<256;i++)
  for (j=0;j<256;j++)
  {
    if (0==abs(random())%50)
      SETBIT(a,i,j);
    else
      CLEARBIT(a,i,j);

    if (0==abs(random())%50)
      SETBIT(b,i,j);
    else
      CLEARBIT(b,i,j);
  }

// Compute binary product for packed matrices
startCPU=CPUtime(); 
packedBinMM(256,a,b,c);
stopCPU=CPUtime();
printf("Row-oriented CPU %.5f\n",stopCPU-startCPU);

for (i=0;i<256;i++)
  for (j=0;j<256;j++)
    if (cSlow[i][j]!=GETBIT(c,i,j))
    {
      printf("c[%d][%d] inconsistency - Row-oriented test failed\n",i,j);
      exit(0);
    }
printf("Row-oriented test passed\n");

// Kronrod's binary matrix multiply
startCPU=CPUtime();
KronrodBinMM(256,a,b,c);
stopCPU=CPUtime();
printf("Kronrod CPU %.5f\n",stopCPU-startCPU);

for (i=0;i<256;i++)
  for (j=0;j<256;j++)
    if (cSlow[i][j]!=GETBIT(c,i,j))
    {
      printf("cSlow[%d][%d]=%d GETBIT=%d inconsistency - Kronrod failed\n",
       i,j,cSlow[i][j],GETBIT(c,i,j));
      exit(0);
    }
printf("Kronrod passed\n");
}