/*
-----------------------------------------------------------------
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 2i of block [*][i] with bit 2i of block 
[*][j], where i is in 0..7, j is in 0..7, and i != j.  Vertical wires
XOR bit 2i+1 of [i][*] with bit 2i+1 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, and maybe complemented.  Hillclimbing was used to find those
constants.  Counting constants are XORed with s to break the symmetry 
between the 64 4x4 blocks, so each differs in at least two ways.

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.

Hillclimbing used 256 seeds created by placing 0..255 in one byte of
one of those 16 64-bit integers, set the rest to all zeros or all 
ones.  All one-bit diffs were tested, and all caused all 1024 bits
to change with non-zero (and non-one) probability after 5 rounds
forward.  That makes these constants unusually good at mixing one-bit
differences with mostly-zero or mostly-one states.  That might make
it unusually bad at something else, who knows.  But my belief is it
makes it better than average at most things.
-----------------------------------------------------------------
*/

#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

#define ROUNDS 256

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

  /* Horizontal linear mixing */
  for (i=0; i<LEN; i+=2) {
    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=1; i<LEN; i+=2) {
    j = x[i];
    j = j ^ (j<<32) ^ (j>>32);
    j = j ^ (j<<16) ^ (j>>48);
    j = j ^ (j<<8) ^ (j>>56);
    x[i] ^= 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[0] ^ y[3] ^ y[7];
  x[ 1] = y[0] ^ y[2] ^ y[4] ^ y[7];
  x[ 2] = y[1] ^ y[4] ^ y[5] ^ y[7];
  x[ 3] = y[0] ^ y[2] ^ y[6] ^ y[7];
  x[ 4] = y[1] ^ y[3] ^ y[4] ^ y[5];
  x[ 5] = y[1] ^ y[2] ^ y[5] ^ y[6];
  x[ 6] = y[0] ^ y[1] ^ y[3] ^ y[4];
  x[ 7] = y[1] ^ y[3] ^ y[6] ^ y[7];
  x[ 8] = (y[ 8]^0x5555555555555555LL) ^ ((y[3] & ~y[5]) | (y[6] & y[0]));
  x[ 9] = (y[ 9]^0x9999999999999999LL) ^ ((y[3] & ~y[4]) | (y[0] & ~y[5]));
  x[10] = (y[10]^0xc3c3c3c3c3c3c3c3LL) ^ ((~y[1] & y[6]) | (y[2] & y[7]));
  x[11] = (y[11]^0xf00ff00ff00ff00fLL) ^ ((~y[3] & y[7]) | (~y[2] & y[1]));
  x[12] = (y[12]^0xff0000ffff0000ffLL) ^ ((~y[4] & ~y[9]) | (~y[3] & ~y[0]));
  x[13] = (y[13]^0xffff00000000ffffLL) ^ ((~y[5] & y[6]) | (y[4] & ~y[1]));
  x[14] = (y[14]^0x00000000ffffffffLL) ^ ((~y[2] & ~y[7]) | (~y[6] & y[5]));
  x[15] = (y[15]^0xffffffffffffffffLL) ^ ((y[4] & ~y[7]) | (~y[1] & ~y[5]));
}

/* 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[2] ^ x[4];
  y[ 1] = x[0] ^ x[1] ^ x[2] ^ x[5] ^ x[7];
  y[ 2] = x[0] ^ x[4] ^ x[5] ^ x[6] ^ x[7];
  y[ 3] = x[0] ^ x[1] ^ x[3] ^ x[6] ^ x[7];
  y[ 4] = x[0] ^ x[3] ^ x[4] ^ x[5];
  y[ 5] = x[0] ^ x[2] ^ x[6];
  y[ 6] = x[0] ^ x[1] ^ x[4] ^ x[5];
  y[ 7] = x[0] ^ x[1] ^ x[2] ^ x[3] ^ x[4] ^ x[6] ^ x[7];
  y[ 8] = (x[ 8]^0x5555555555555555LL) ^ ((y[3] & ~y[5]) | (y[6] & y[0]));
  y[ 9] = (x[ 9]^0x9999999999999999LL) ^ ((y[3] & ~y[4]) | (y[0] & ~y[5]));
  y[10] = (x[10]^0xc3c3c3c3c3c3c3c3LL) ^ ((~y[1] & y[6]) | (y[2] & y[7]));
  y[11] = (x[11]^0xf00ff00ff00ff00fLL) ^ ((~y[3] & y[7]) | (~y[2] & y[1]));
  y[12] = (x[12]^0xff0000ffff0000ffLL) ^ ((~y[4] & ~y[9]) | (~y[3] & ~y[0]));
  y[13] = (x[13]^0xffff00000000ffffLL) ^ ((~y[5] & y[6]) | (y[4] & ~y[1]));
  y[14] = (x[14]^0x00000000ffffffffLL) ^ ((~y[2] & ~y[7]) | (~y[6] & y[5]));
  y[15] = (x[15]^0xffffffffffffffffLL) ^ ((y[4] & ~y[7]) | (~y[1] & ~y[5]));

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

  /* Horizontal linear mixing */
  for (i=0; i<LEN; i+=2) {
    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=1; i<LEN; i+=2) { 
    j = x[i];
    j = j ^ (j<<32) ^ (j>>32);
    j = j ^ (j<<16) ^ (j>>48);
    j = j ^ (j<<8) ^ (j>>56);
    x[i] ^= 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");
}

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;
}

void driver(ranctx *rctx)
{
  u8 x[LEN], y[LEN];
  u4 m,n,o,p,count,i,j,k;
  float val, worst = 1.0;
  float zcount = 0.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) {
	for (p=0; p<2*LEN; ++p) {
	  float c[ROWS][COLS][LEN];  /* count of how many times bit flipped */
	  u4 r;  /* count # of rounds */

	  /* if (m != 0 || n != 0 || o != 0 || p != 0) continue; */
	  /* if (m != 0) continue; */

	  /* 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); */
	      x[i] = y[i] = (p/LEN ? ~((u8)0) : ((u8)0));
	    }
	    x[p%LEN] ^= ((u8)r)<<0;
	    y[p%LEN] = x[p%LEN];

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

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

	    ihash(x); ihash(y);
	    ihash(x); ihash(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 */
	  val = showc(c, 0);
	  if (val == 0) {
	    ++zcount;
	  }
	  if (val < worst) {
	    worst = val;
	    printf("new worst: %f  %d %d %d\n", worst, m,n,o);
	  }
	}
      }
    }
  }
}

int main()
{
  int i;
  ranctx rctx = {1,1,1,63};
  for (i=0; i<20; ++i) ranval(&rctx);
  driver(&rctx);
}