/*
-----------------------------------------------------------------
Software simulation of a 1024-bit hardware hash.  A work-in-progress
with the goal of entering the NIST cryptographic hash competition.

Bits are arranged in a 32 x 32 grid, well, an 8x8 grid of 
4x4-bit blocks.  The hash consists of horizontal and vertical
XORs, followed by the a nonlinear mix applied to every 4x4-bit block.

The horizontal wires XOR bit i of block [*][i] with bit i of block 
[*][j], where i in 0..7, j in 0..7, and i != j.  Vertical wires
XOR bit i+8 of [i][*] with bit i+8 of block [j][*] where i in 0..7,
j in 0..7, and i != j.

The nonlinear mix of 4x4 blocks is actually an LFSR for bits 0..7,
and (!s)^((w&x)|(y&z)) for bits 8..15, where w,x,y,z are all less 
than s.

The hash is done in rounds.  Each round does the horizontal and 
vertical XORs, then mixes the 4x4 blocks.

The internal state is represented by 16 64-bit integers, x[16],
where x[k] & (1<<(8*i+j)) is the kth bit of block [i][j].  This 
allows the nonlinear mixing to be done for all blocks simulateously,
and also allows horizontal and vertical wires to be handled in bulk.
-----------------------------------------------------------------
*/

#include <stdio.h>
#include <stddef.h>

typedef unsigned long  u4;
typedef unsigned long long u8;

typedef struct ranctx { u4 a; u4 b; u4 c; u4 d; } ranctx;

#define rot(x,k) ((x<<k)|(x>>(32-k)))
u4 ranval( ranctx *x ) {
  u4 e = x->a - rot(x->b, 27);
  x->a = x->b ^ rot(x->c, 17);
  x->b = x->c + x->d;
  x->c = x->d + e;
  x->d = e + x->a;
  return e;
}

#define LEN 16

#define ROWS 8
#define COLS 8


/* do one round of the hash */
void hash(u8 x[LEN])
{
  u8 i, j;
  u8 y[LEN];

  /* Horizontal linear mixing */
  for (i=0; i<8; ++i) {
    j = x[i];
    j = j ^ ((j&0x5555555555555555LL)<<1) ^ ((j&0xaaaaaaaaaaaaaaaaLL)>>1);
    j = j ^ ((j&0x3333333333333333LL)<<2) ^ ((j&0xccccccccccccccccLL)>>2);
    j = j ^ ((j&0x0f0f0f0f0f0f0f0fLL)<<4) ^ ((j&0xf0f0f0f0f0f0f0f0LL)>>4);
    x[i] ^= j;
  }

  /* Vertical linear mixing */
  for (i=0; i<8; ++i) {
    j = x[i+8];
    j = j ^ (j<<32) ^ (j>>32);
    j = j ^ (j<<16) ^ (j>>48);
    j = j ^ (j<<8) ^ (j>>56);
    x[i+8] ^= j;
  }

  for (i=0; i<LEN; ++i) y[i] = x[15-i];

  /* Do the nonlinear mixing of the bits in the 4x4 blocks */
  x[ 0] = y[2] ^ y[3] ^ y[6] ^ y[7];
  x[ 1] = y[0] ^ y[1] ^ y[2] ^ y[3];
  x[ 2] = y[0] ^ y[5] ^ y[6] ^ y[7];
  x[ 3] = y[0] ^ y[1] ^ y[2] ^ y[5];
  x[ 4] = y[1] ^ y[2] ^ y[3] ^ y[4];
  x[ 5] = y[0] ^ y[4] ^ y[6] ^ y[7];
  x[ 6] = y[0] ^ y[1] ^ y[2] ^ y[6];
  x[ 7] = y[0] ^ y[2] ^ y[3];
  x[ 8] = (~y[ 8]) ^ ((y[7] & y[6]) | (y[3] & y[0]));
  x[ 9] = (~y[ 9]) ^ ((y[8] & y[6]) | (y[7] & y[2]));
  x[10] = (~y[10]) ^ ((y[6] & y[1]) | (y[4] & y[3]));
  x[11] = (~y[11]) ^ ((y[8] & y[3]) | (y[4] & y[1]));
  x[12] = (~y[12]) ^ ((y[11] & y[2]) | (y[10] & y[1]));
  x[13] = (~y[13]) ^ ((y[12] & y[0]) | (y[6] & y[5]));
  x[14] = (~y[14]) ^ ((y[12] & y[2]) | (y[10] & y[4]));
  x[15] = (~y[15]) ^ ((y[11] & y[5]) | (y[7] & y[4]));
}


/* do the inverse of one round of the hash */
void ihash(u8 x[LEN])
{
  u8 i, j;
  u8 y[LEN];

  /* Do the nonlinear mixing of the bits in the 4x4 blocks */
  y[ 0] = x[0] ^ x[1] ^ x[4] ^ x[5] ^ x[7];
  y[ 1] = x[1] ^ x[7];
  y[ 2] = x[2] ^ x[3] ^ x[4] ^ x[5] ^ x[7];
  y[ 3] = x[0] ^ x[1] ^ x[2] ^ x[3] ^ x[7];
  y[ 4] = x[0] ^ x[5] ^ x[7];
  y[ 5] = x[0] ^ x[2] ^ x[7];
  y[ 6] = x[0] ^ x[2] ^ x[3] ^ x[6] ^ x[7];
  y[ 7] = x[0] ^ x[1] ^ x[2] ^ x[3] ^ x[4] ^ x[5] ^ x[6] ^ x[7];
  y[ 8] = (~x[ 8]) ^ ((y[7] & y[6]) | (y[3] & y[0]));
  y[ 9] = (~x[ 9]) ^ ((y[8] & y[6]) | (y[7] & y[2]));
  y[10] = (~x[10]) ^ ((y[6] & y[1]) | (y[4] & y[3]));
  y[11] = (~x[11]) ^ ((y[8] & y[3]) | (y[4] & y[1]));
  y[12] = (~x[12]) ^ ((y[11] & y[2]) | (y[10] & y[1]));
  y[13] = (~x[13]) ^ ((y[12] & y[0]) | (y[6] & y[5]));
  y[14] = (~x[14]) ^ ((y[12] & y[2]) | (y[10] & y[4]));
  y[15] = (~x[15]) ^ ((y[11] & y[5]) | (y[7] & y[4]));

  for (i=0; i<LEN; ++i) x[15-i] = y[i];

  /* Horizontal linear mixing */
  for (i=0; i<8; ++i) {
    j = x[i];
    j = j ^ ((j&0x5555555555555555LL)<<1) ^ ((j&0xaaaaaaaaaaaaaaaaLL)>>1);
    j = j ^ ((j&0x3333333333333333LL)<<2) ^ ((j&0xccccccccccccccccLL)>>2);
    j = j ^ ((j&0x0f0f0f0f0f0f0f0fLL)<<4) ^ ((j&0xf0f0f0f0f0f0f0f0LL)>>4);
    x[i] ^= j;
  }

  /* Vertical linear mixing */
  for (i=0; i<8; ++i) {
    j = x[i+8];
    j = j ^ (j<<32) ^ (j>>32);
    j = j ^ (j<<16) ^ (j>>48);
    j = j ^ (j<<8) ^ (j>>56);
    x[i+8] ^= j;
  }
}




void show(u8 x[LEN])
{
  int i,j,k;
  for (i=0; i<ROWS; ++i) {
    for (j=0; j<COLS; ++j) {
      for (k=0; k<4; ++k) {
	printf(" %s", (x[k] & (((u8)1)<<(COLS*i+j))) ? "1" : "0");
      }
    }
    printf("\n");
    for (j=0; j<COLS; ++j) {
      for (k=8; --k>=4; ) {
	printf(" %s", (x[k] & (((u8)1)<<(COLS*i+j))) ? "1" : "0");
      }
    }
    printf("\n");
    for (j=0; j<COLS; ++j) {
      for (k=8; k<12; ++k) {
	printf(" %s", (x[k] & (((u8)1)<<(COLS*i+j))) ? "1" : "0");
      }
    }
    printf("\n");
    for (j=0; j<COLS; ++j) {
      for (k=16; --k>=12; ) {
	printf(" %s", (x[k] & (((u8)1)<<(COLS*i+j))) ? "1" : "0");
      }
    }
    printf("\n");
  }
  printf("\n");
}

#define ROUNDS 10000

void showx(float x, float *worst, int verbose)
{
  x = x/ROUNDS;
  if (x < *worst) *worst = x;
  if (verbose) printf("%3d", (int)(100*x + 0.5));
  if (1.0-x < *worst) *worst = 1.0-x;
}

float showc(float c[ROWS][COLS][LEN], int verbose)
{
  int i,j,k;
  float worst = 1.0;
  for (i=0; i<ROWS; ++i) {
    for (j=0; j<COLS; ++j) {
      for (k=0; k<4; ++k) {
	showx(c[i][j][k], &worst, verbose);
      }
      if (verbose) printf(" ");
    }
    if (verbose) printf("\n");
    for (j=0; j<COLS; ++j) {
      for (k=8; --k>=4; ) {
	showx(c[i][j][k], &worst, verbose);
      }
      if (verbose) printf(" ");
    }
    if (verbose) printf("\n");
    for (j=0; j<COLS; ++j) {
      for (k=8; k<12; ++k) {
	showx(c[i][j][k], &worst, verbose);
      }
      if (verbose) printf(" ");
    }
    if (verbose) printf("\n");
    for (j=0; j<COLS; ++j) {
      for (k=16; --k>=12; ) {
	showx(c[i][j][k], &worst, verbose);
      }
      if (verbose) printf(" ");
    }
    if (verbose) printf("\n\n");
  }
  /* printf("  worst = %f\n", worst); */
  return worst;
}

int main()
{
  ranctx rctx = {1,1,1,1};
  u8 x[LEN], y[LEN];
  u4 m,n,o,count,i,j,k;
  float val, worst = 1.0;

  /* see how thoroughly bit "tab[m][n] & (1<<o)" propogates */
  for (m=0; m<ROWS; ++m) {
    for (n=0; n<COLS; ++n) {
      for (o=0; o<LEN; ++o) {
	/* if (m != 6 || n != 1 || o != 11) continue; */

	float c[ROWS][COLS][LEN];  /* count of how many times bit flipped */
	u4 r;  /* count # of rounds */

	/* initialize c */
	for (i=0; i<ROWS; ++i) {
	  for (j=0; j<COLS; ++j) {
	    for (k=0; k<LEN; ++k) {
	      c[i][j][k] = 0;
	    }
	  }
	}

	/* try r keys */
	for (r=0; r<ROUNDS; ++r) {

	  /* initialize to a random key */
	  for (i=0; i<LEN; ++i) {
	    x[i] = y[i] = ((u8)ranval(&rctx)) | (((u8)ranval(&rctx))<<32);
	  }

	  /* start with bit [m][n][o] different */
	  y[o] ^= (((u8)1) << (COLS*m+n));

	  /* apply several rounds of the hash to both keys */
	  hash(x); hash(y);
	  hash(x); hash(y);
	  hash(x); hash(y);

	  hash(x); hash(y);
	  hash(x); hash(y);

	  /* measure which bits were affected */
	  for (k=0; k<LEN; ++k) {
	    u8 diff = x[k] ^ y[k];
	    for (i=0; i<ROWS; ++i) {
	      for (j=0; j<COLS; ++j) {
		if (diff & (((u8)1) << (COLS*i+j))) {
		  ++c[i][j][k];
		}
	      }
	    }
	  }
	}

	/* report which bits were affected */
	/* printf("%2d %2d %2d:\n", m,n,o); */
	val = showc(c, 0);
	if (val < worst) {
	  worst = val;
	  printf("new worst: %f  %d %d %d\n", worst, m,n,o);
	}
      }
    }
  }
}