rank.c

随机数算法

/*
	Compute rank of square matrix from raw bit data.
	Enter with source pointer into data, dimension of matrix and endianness.
	Returns with rank of matrix and source pointer updated to next batch of
	bits.  Returns -1 if requested dimension is larger than MAXRANK.

	Definition of rank: "The maximum number of linearly independent row or
	column vectors.  Theorem 14: For any matrix A, the determinant rank, 
	row rank and column rank are equal." pg 511, "Advanced Engineering
	Mathmatics", C. Ray Wylie, 1975

	Author = Mike Rosing
	 Date = Nov. 2, 2001
*/

#include "randtest.h"

#define dimchar	(MAXRANK/8)

int binrank( SRCPTR  *src, int dim, int endian)
{
	int row, rankval, j, k, i, dumdim, mscount;
	unsigned char matrix[MAXRANK][dimchar+1];
	unsigned char bitmask, save;
	int ji, ki, ii, rowi;
	unsigned char bitmaski;
	
	if( dim > MAXRANK) return (-1);

/*  clear out local storage  */

	for(row=0; row<MAXRANK; row++)
	  for(j=0; j<dimchar; j++)
	    matrix[row][j] = 0;

/*  copy data from source to matrix format.
	this makes bit twiddling easier.
*/

/*	for(i=0; i<=(dim*dim)/8; i++)
	{
	  printf("%x ",src->byteptr[i]);
	  if (!((i&3) ^ 3)) printf("\n");
	} 
	printf("\n");
//	printf("starting address = %x\n", src->byteptr);
*/
	for( row=0; row<dim; row++) getbits( src, matrix[row], endian, dim);
	mscount = dim/8;
	if( dim&7) mscount++;

/*  print out matrix data in ascii format for human eye to deal with */

	dumdim = dim - 1;

/*	for(row=0; row<dim; row++)
	{
	  bitmask = 1 << ((dumdim) & 7);
	  k = 0;
	  for(j=0; j<dim; j++) 
	  {
	    if( !endian)  k = (dumdim-j)/8;
//	    printf("j=%d k=%d bitmask= %x\n", j,k, bitmask);
	    i = bitmask & matrix[row][k] ? 1 : 0;
	    printf("%x ", i);
	    bitmask >>= 1;
	    if (!bitmask)
	    { 
	      bitmask = 0x80;
	      if(endian) k++;
	    }
	  }
	  printf("\n");
//	  for(i=0; i<=dim/8; i++) printf("%x ", matrix[row][i] );
//	  printf("\n");
	}
*/
/*  clear out matrix.  Number of 1's on diagonal tells us rank
	(same as the number of independent vectors in the matrix
         but we aren't solving anything, just counting)  */
	
	for( row=0; row<dim; row++)
	{
	
/*  set up byte and bit to be on diagonal  */

	        dumdim = dim - 1 - row;
		bitmask = 1 << ( dumdim & 7);
		j = dumdim/8;
		if( endian)
		{
	          if( dim<8 ) j=0;
		  else if( row < (dim&7)) j = 0;
		  else 
		  {
		    j = (row - (dim&7))/8 ;
		    if( dim&7) j++;
		  }
		}

/*  look for bit set on diagonal.  Swap rows if found */

//		printf("looking for diagonal, j = %d bitmask = %x\n", j, bitmask);
		if( !(bitmask & matrix[row][j]))
		{
			for( k=row+1; k<dim; k++)
			{
				if( bitmask & matrix[k][j])
				{
//				  printf("swapping row %d with %d\n", row, k);
					for( i=0; i<=dimchar; i++)
					{
						save = matrix[k][i];
						matrix[k][i] = matrix[row][i];
						matrix[row][i] = save;
					}
					break;
				}
			}
		}

/*  double check.  If found, clear bits below this row
	otherwise skip to next row (column is empty)
*/

		if( bitmask & matrix[row][j])
		{
			for( k=row+1; k<dim; k++)
			{
			  if( bitmask & matrix[k][j])
			    {
			      for( i=0; i<=dim/8; i++)
				matrix[k][i] ^= matrix[row][i];
			    }
			}
		}
	}// end of row loop
/*printf("\n");

	for(row=0; row<dim; row++)
	{
	  bitmask = 1 << ((dim-1) & 7);
	  k = 0;
	  for(j=0; j<dim; j++) 
	  {
	    if( !endian)  k = (dim-1-j)/8;
	    i =  bitmask & matrix[row][k] ? 1 : 0;
	    printf("%x ",i);
	    bitmask >>= 1;
	    if (!bitmask)
	    {
	      bitmask = 0x80;
	      if(endian) k++;
	    }
	  }
	  printf("\n");
	}
*/
/*  next do back solve portion
    The main thing we need to do is remove rows that are 
    duplicates.
*/

	for( row=dim-1; row>=0; row--)
	{
	
/*  set up byte and bit to be on diagonal  */

	        dumdim = dim - 1 - row;
		bitmask = 1 << ( dumdim & 7);
		j = dumdim/8;
		if( endian)
		{
	          if( dim<8 ) j=0;
		  else if( row < (dim&7)) j = 0;
		  else 
		  {
		    j = (row - (dim&7))/8 ;
		    if( dim&7) j++;
		  }
		}

/*  look for bit set on diagonal.  Swap rows if found */

		if( !(bitmask & matrix[row][j]))
		{
			for( k=row-1; k>=0; k--)
			{
				if( bitmask & matrix[k][j])
				{
					for( i=0; i<=dimchar; i++)
					{
						save = matrix[k][i];
						matrix[k][i] = matrix[row][i];
						matrix[row][i] = save;
					}
					break;
				}
			}
		}

/*  double check.  If found, clear bits above this row
	otherwise skip to next row (column is empty)
*/

		if( bitmask & matrix[row][j])
		{
			for( k=row-1; k>=0; k--)
			{
			  if( bitmask & matrix[k][j])
			    {
			      for( i=0; i<=dim/8; i++)
				matrix[k][i] ^= matrix[row][i];
			    }
			}
		}
	}// end of row loop

/*  count the number of non zero rows  */

	rankval = 0;
	for( row=0; row<dim; row++)
	{
	  for (j=0; j<mscount; j++)
	  { 
	    if (matrix[row][j])
	    {
	      rankval++;
	      break;
	    }
	  }
	}

//        printf("\n");

/*  print out matrix data in ascii format for human eye to deal with */

/*
	for(row=0; row<dim; row++)
	{
	  bitmask = 1 << ((dim-1) & 7);
	  k = 0;
	  for(j=0; j<dim; j++) 
	  {
	    if( !endian)  k = (dim-1-j)/8;
	    i =  bitmask & matrix[row][k] ? 1 : 0;
	    printf("%x ",i);
	    bitmask >>= 1;
	    if (!bitmask)
	    {
	      bitmask = 0x80;
	      if(endian) k++;
	    }
	  }
	  printf("\n");
	}
*/
	return rankval;
}


/*  compute 1/2^x for genMatrixP routine below */

double invPow( int x)
{
  double num;
  int i;

  if( !x) return 1.0;

  num = 1.0;
  for( i=0; i<x; i++) num /= 2.0;
  return num;
}

/*  generate probability matrix for any number of square
    matricies.  Formula from NIST page 70 and its references:

                               r-1
    p(r,M) = 1/2^((M-r)^2) * product{ (1 - 2^(i-M))^2 /(1 - 2^(i-r))}
                               i=0

*/
/*
void genMatrixP( int maxdim)
{
  double pout, ptop, pbot, sum;
  int M, r, i;

  printf("\n M  r  probability\n\n");
  for( M=1; M<=maxdim; M++)
  {
    sum = 0.0;
    for( r=0; r<=M; r++)
    {
      pout = invPow( (M-r)*(M-r));
      for( i=0; i<r; i++)
      {
	ptop = 1.0 - invPow(M-i);
	pbot = 1.0 - invPow(r-i);
	pout *= ptop*ptop/pbot;
      }
      sum += pout;
      printf("%3d %3d %e\n", M, r, pout);
    }
    printf("sum = %f\n\n", sum);
  }
}
*/

/*  Same thing, but one at a time.  */

double genMatrixP( int M, int r)
{
  double pout, ptop, pbot;
  int i;

  pout = invPow( (M-r)*(M-r));
  for( i=0; i<r; i++)
  {
    ptop = 1.0 - invPow(M-i);
    pbot = 1.0 - invPow(r-i);
    pout *= ptop*ptop/pbot;
  }
  return pout;
}

/*  Binary matrix rank test.

    Described in NIST and used in DIEHARD the following test is designed for
    megabits of data.  It is assumed that the matrix dimension is > 8
    for sensible results.  The math should be correct for any dimension,
    but this code is hard wired to deal with dimesion >= 6 so the actual
    cut off is 8.  The cardinality of 5 ranks is computed and the chi^2
    observed is then found using exact formulas for the probabilities.
    The p-value is found from 1 - P(chi^2, 5) since we have 5 degrees
    of freedom.

    Enter with source pointer, endianness, matrix dimension and total
      number of matricies to deal with.
    Returns chi^2(observed) and p-value, function is negative on error,
      positive on success (negative function return means chi^2 and
      p-value should be ignored).

*/

int ranktest( SRCPTR src, int endian, int Ndim, int Mtodo,
	      double *chsqr, double *pvalue)
{
  double F[6], P[6], t, sum;

 /*  The first 5 elements in each array correspond to ranks of
     M, M-1, M-2, M-3 and M-4 respectively.  The last element
     contains the sum of all remaining ranks (and should be zero
     for truely random data on this scale of less than 1e6 matricies).
 */

  int i, j, k;
  SRCPTR bits;

  if( Ndim < 8) return -1;
  bits.byteptr = src.byteptr;
  bits.bitpos = src.bitpos;

/*  initialize counts to zero  */

  for( i=0; i<6; i++) F[i] = 0;

/*  find rank of matrix and increment correct bin  */

  for( i=0; i<Mtodo; i++)
  {
    k = binrank( &bits, Ndim, endian);
    if( k<0)
    {
      printf("error computing binrank for dimension %d\n", Ndim);
      return -1;
    }
    if( k> Ndim-5) F[Ndim-k]++;
    else F[5]++;
  }

/*  compute probabilities for each rank  */

  sum = 0.0;
  for( i=0; i<5; i++)
  {
    P[i] = genMatrixP( Ndim, Ndim-i);
    sum += P[i];
//    printf("P[%d] = %f  F[%d] = %f\n", i, P[i], i, F[i]);
  }
  P[5] = 1.0 - sum;
//  printf("P[5] = %e  F[5] = %f\n", P[5], F[5]);

/*  compute chi^2(observed) and p-value  */

  sum = 0.0;
  j = 0;
  for( i=0; i<5; i++)
  {
    j += F[i];
    t = F[i] - P[i]*Mtodo;
//    printf("i= %d  t = %f\n", i, t);
    sum += t*t/P[i]/Mtodo;
  }
  t = (Mtodo - j) - P[5]*Mtodo;
  sum += t*t/P[5]/Ndim;
  *pvalue = chisqr( sum, 5);
//  printf("chisqr = %f, pvalue = %f \n", sum, *pvalue);
  *chsqr = sum;
  return 1;
}