randtest.c

随机数算法

/* randtest.c : Defines the entry point for the console application.
//
//
//	RANDTEST is a command-line program.  It is designed so that it can be executed
//	from a batch file.  It requires the name of an input file and the word size
//	(in bits) to be used during the analysis.  If the input file is supplied on the
//	command line, it must be the first argument.  If not, it will be requested.
//	The word size may be specified on the command line using the -wNN option.
//	If it is not specified, it will be requested.
/*
	This file originally written by John Hadstate  
	Modified by Mike Rosing on and after October 19, 2001
		added fixed ram for data block and a subroutine to 
		access it as binary.
		added -x option for binary input file (for rank test)

		added -r option for binary matrix rank test (Nov 10, 2001)
		added -m option for monobits frequency test (Nov 14, 2001)
                added -l option for longest run of ones test (Jan 2, 2002)
                      This includes bigfloat package for computing probability
		      tables to high accuracy for low runs.
*/
#include "randtest.h"
#include "bigfloat.h"

#define MAXRAMSIZE 16*1024*1024

int nByteMask[BITS_PER_BYTE] =
{
	0x00000001, 0x00000002, 0x00000004, 0x00000008,
	0x00000010, 0x00000020, 0x00000040, 0x00000080
};

unsigned char rawbits[MAXRAMSIZE];

int main(int argc, char* argv[])
{
	FILE		*fpInput = NULL;

	int			nBitNumber = (-1),
				nWordSize = (-1);
	int 		i;
	char		*pc = NULL;
	char		strOptions[1024];
	char		strInputFile[1024];

	bool		bAnalyzeAllBits = false;
	bool		bNumericAnalysis = false;
	bool            dorank = false;
	bool            domono = false;
	bool            dofreq = false;
	bool            doruns = false;
	bool            dolongs = false;

	int		nBytes;
	SRCPTR	        bitBlockltl, bitBlockbig;
	int		rankvalue, ii, rankstart;
	int             error;
	double          chi, pvalue;
	int             monowidth;
	int             numPs, numBlocks;
	int             longRunBlocks, tblsz;

	CHIENTRY       ctbl[21];
	FLOAT          test;

/*  initialize floating point package  */

	init_float();

/*	BuildChiTable( 8, ctbl);
	printf("M = 8\n");
	for( i=0; i<8; i++)
	  printf("ctbl[%d].class = %d  .prob = %e\n",i, ctbl[i].class, ctbl[i].prob);
	printf("\n");

	printf("M = 128\n");
	BuildChiTable( 128, ctbl);
	i = 0;
	for( i=0; i<21; i++)
	  printf("ctbl[%d].class = %d  .prob = %e\n",i, ctbl[i].class, ctbl[i].prob);
	printf("\n");
/*
	BuildChiTable( 512, ctbl);
	for( i=0; i<21; i++)
	  printf("ctbl[%d].class = %d  .prob = %e\n",i, ctbl[i].class, ctbl[i].prob);

	exit(0);       
/*
	First, concatenate all the command-line options that are
	prefixed with a '-' and convert them to lower case.
	We are pretty sure these cannot be the Input File name
	and, in any case, the Input File name may not start with
	a '-' (our restriction).
*/
	strOptions[0] = NUL;

	for (i = 1; i < argc; i++)
	{
		if (*argv[i] == '-') SAFECAT(strOptions, argv[i]);
	}

	for (i = 0; strOptions[i] != 0; i++) strOptions[i] = tolower(strOptions[i]);
/*
	Now go find the first command-line option that does not begin with
	a '-'.  Accept this as the Input File name.
*/
	strInputFile[0] = NUL;

	for (i = 1; i < argc; i++)
	{
		if (*argv[i] != '-')
		{
//			SAFECOPY(strInputFile, argv[1]);  I think this is a bug!! (MGR)
			SAFECOPY(strInputFile, argv[i]);
			break;
		}
	}
/*
	If the Input File name has not been specified, we've got to
	prompt for it interactively.  We're going to read from stdin
	to allow one last bit of hacking: stdin could be redirected to
	a file containing the Input File name and other required
	parameters.
*/

	if (strInputFile[0] == NUL)
	{
		fputs("\nName of input file: ", stdout);
		fscanf(stdin, " %1023s", strInputFile);
	}
/*
	If the word size of the random number stream has not been specified,
	we've got to prompt for it interactively.
*/
	if ((pc = strstr(strOptions, "-w")) == NULL)
	{
		fputs("\nWord size (in bits): ", stdout);
		fscanf(stdin, " %d", &nWordSize);
	}
	else
	{
		if (sscanf(pc, " -w%d", &nWordSize) != 1)
		{
			fprintf(stdout, "\n-w Option Specification Invalid: %s\n", pc);

			return (1);
		}
	}
/*
	If there are no command line options specified,
	we need to prompt for numeric or bitwise analysis.
*/
	if (strOptions[0] == NUL)
	{
		char		cTemp;

		fputs("\nNumeric (n) or Bitwise (b) analysis? ", stdout);
		fscanf(stdin, " %c", &cTemp);
		cTemp = tolower(cTemp);

		if (cTemp == 'n')
		{
			SAFECOPY(strOptions, "-n");
		}
		else if (cTemp == 'b')
		{

			fputs("\nWhich bit (-1 for ALL): ", stdout);
			fscanf(stdin, " %d", &nBitNumber);

			if (nBitNumber >= 0)
			{
				sprintf(strOptions, "-b%d", nBitNumber);
			}
		}
	}
/*
	If the -n option is specified, we don't need the -b option.  If the
	-b option is needed but not specified, it means run the analysis once
	on each bit in the word.
*/
	if ((pc = strstr(strOptions, "-n")) != NULL)
	{
		if ((pc = strstr(strOptions, "-b")) != NULL)
		{
			fprintf(stdout, "\n-n Option conflicts with -b Option: %s\n", strOptions);

			return (1);
		}

		bNumericAnalysis = true;
	}
	else
	{
		bNumericAnalysis = false;
	/*
		If the bit number in the words of the random number stream has not been specified,
		we assume that they want to analyze all bits.
	*/
		if ((pc = strstr(strOptions, "-b")) != NULL)
		{
			if (sscanf(pc, " -b%d", &nBitNumber) != 1)
			{
				fprintf(stdout, "\n-b Option Specification Invalid: %s\n", pc);

				return (1);
			}

			bAnalyzeAllBits = false;
		}
		else
		{
			bAnalyzeAllBits = true;
		}
	}
/*
	Word Size must be greater than 0 bits.  For numeric analysis, word size
	must be less than or equal to 32 bits.
*/
	if ((nWordSize <= 0) || (bNumericAnalysis && (32 < nWordSize)))
	{
		fprintf(stdout, "\n-w Option Size Invalid: %d\n", nWordSize);

		return (1);
	}
/*
	If Bit Number is required, it must be between 0 and (N-1) where N
	is the number of bits in a word.
*/
	if (!bNumericAnalysis && !bAnalyzeAllBits)
	{
		if ((nBitNumber < 0) || (nWordSize <= nBitNumber))
		{
			fprintf(stdout, "\n-b Option Bit Number Invalid: %d\n", nBitNumber);

			return (1);
		}
	}

	if ((fpInput = fopen(strInputFile, "r")) == NULL)
	{
		char	strErrorMessage[1024];

		sprintf(strErrorMessage, "\nUnable to open %s", strInputFile);
		perror(strErrorMessage);

		return (2);
	}

/*  Read up to max buffer size of data.  We can loop on this if file is
	bigger, or increase buffer size.  Ram is cheap nowadays!!
*/
	if( ( pc = strstr( strOptions, "-x") ) != NULL)
	{
		nBytes = fread( rawbits, 1, MAXRAMSIZE, fpInput);
		bitBlockltl.byteptr = rawbits;
		bitBlockltl.bitpos = 0;  // little endian, least significant bit
		bitBlockbig.byteptr = rawbits;
		bitBlockbig.bitpos = 7;  // big endian, most significant bit
	}

/*  if rank test specified, get starting matrix dimension.
    Test loops over data block going from starting dimension up to
    maximum word size.
*/
	if( (pc = strstr( strOptions, "-r") ) != NULL)
	{
	  if (sscanf(pc, " -r%d", &rankstart) != 1) 
	  {
	    fprintf(stdout, "\ninvalid parameter for rank matrix probabilities\n");
	    return(3);
	  }
	  dorank = true;
	}

/*  if monobits test specified, get number of bits in block to test.
    does as many blocks as possible with raw data.
*/
	if( (pc = strstr( strOptions, "-m") ) != NULL)
	{
	  if (sscanf(pc, " -m%d", &monowidth) != 1) 
	  {
	    fprintf(stdout, "\ninvalid parameter for monobits frequency test\n");
	    return(4);
	  }
	  domono = true;
	}

/*  if frequency test specified, get number of p-values to generate.
    does as many blocks as possible with raw data, assuming that nWordSize
    is bits per block.  Number of bits is assumed to be 8*nBytes, so number
    of bits to work with is 8*nBytes/(number of p-values)(=N in NIST docs)
    and number of trials within that is determined by nWordSize (=M in 
    NIST docs).
*/
	if( (pc = strstr( strOptions, "-f") ) != NULL)
	{
	  if (sscanf(pc, " -f%d", &numPs) != 1) 
	  {
	    fprintf(stdout, "\ninvalid parameter for frequency test within a block\n");
	    return(4);
	  }
	  dofreq = true;
	}

/*  if runs test specified, get number of p-values to generate.
    Does as many blocks as possible assuming number of bits is 8*nBytes.
*/
	if( (pc = strstr( strOptions, "-u") ) != NULL)
	{
	  if (sscanf(pc, " -u%d", &numPs) != 1) 
	  {
	    fprintf(stdout, "\ninvalid parameter for frequency test within a block\n");
	    return(4);
	  }
	  doruns = true;
	}

/*  if longest runs test specified, get number of p-values to generate.
    Does as many blocks as possible assuming number of bits is 8*nBytes.
*/
	if( (pc = strstr( strOptions, "-l") ) != NULL)
	{
	  if (sscanf(pc, " -l%d", &longRunBlocks) != 1) 
	  {
	    fprintf(stdout, "\ninvalid parameter for frequency test within a block\n");
	    return(5);
	  }
	  dolongs = true;
	}

// need to merge both methods.  Rank test assumes binary data already
// internal.  compare running times and simplicity, but merge both methods.		

	if (bNumericAnalysis)
	{
	//	AnalyzeNumeric(fpInput, nWordSize);
	}
	else
	{
	  if (bAnalyzeAllBits)
	  {
	    if (dorank)
	    {
	      for (ii = rankstart; ii <= nWordSize; ii++)
	      {
		
/*  test matrix rank code  */

		printf("\ndimension = %d\n", ii);
		error = ranktest( bitBlockltl, LITTLENDIAN, ii, 8*nBytes/ii/ii,
				  &chi, &pvalue);
		if( error < 0) printf("Error computing little endian binary rank test.\n");
		else printf("Little endian chi^2 = %f,  p-value = %f\n", chi, pvalue);
		error = ranktest( bitBlockbig, BIGENDIAN, ii, 8*nBytes/ii/ii,
				  &chi, &pvalue);
		if( error < 0) printf("Error computing big endian binary rank test.\n");
		else printf("Big endian chi^2 = %f,  p-value = %f\n", chi, pvalue);
	      }
	    }
	    if (domono)
	    {
	      ii = nBytes*8/monowidth;
	      for( i=0; i<ii; i++)
	      {
		chi = monobits( &bitBlockltl, monowidth);
		printf("block %d gives p-value %f\n", i, chi);
	      }
	    }
	    if (dofreq)
	    {
      /* reset source pointers */

	      bitBlockltl.byteptr = rawbits;
	      bitBlockltl.bitpos = 0;  // little endian, least significant bit
	      bitBlockbig.byteptr = rawbits;
	      bitBlockbig.bitpos = 7;  // big endian, most significant bit
	      numBlocks = 8*nBytes/numPs/nWordSize;
	      printf("numblocks = %d\n", numBlocks);
	      for( i=0; i<numPs; i++)
	      {
		chi = freqtest( &bitBlockltl, LITTLENDIAN, nWordSize, numBlocks);
		pvalue = chisqr( chi, numBlocks);
		printf("frequency test little endian chi^2 = %f, p-value = %f\n", 
		       chi, pvalue);
		chi = freqtest( &bitBlockbig, BIGENDIAN, nWordSize, numBlocks);
		pvalue = chisqr( chi, numBlocks);
		printf("frequency test big endian chi^2 = %f, p-value = %f\n", 
		       chi, pvalue);
	      }
	    }
	    if( doruns)
	    {
	      bitBlockltl.byteptr = rawbits;
	      bitBlockltl.bitpos = 0;  // little endian, least significant bit
	      bitBlockbig.byteptr = rawbits;
	      bitBlockbig.bitpos = 7;  // big endian, most significant bit
	      numBlocks = 8*nBytes/numPs;
	      printf("number of bits = %d\n", numBlocks);
	      for( i=0; i<numPs; i++)
	      {
		pvalue = runs( &bitBlockltl, LITTLENDIAN, numBlocks );
		printf("runs test little endian p-value = %f\n",  pvalue);
		pvalue = runs( &bitBlockbig, BIGENDIAN, numBlocks);
		printf("runs test big endian p-value = %f\n", pvalue);
	      }
	    }
	    if( dolongs)
	    {
	      bitBlockltl.byteptr = rawbits;
	      bitBlockltl.bitpos = 0;  // little endian, least significant bit
	      bitBlockbig.byteptr = rawbits;
	      bitBlockbig.bitpos = 7;  // big endian, most significant bit
	      printf("WARNING: large block sizes may take a very long time to simply\n");
	      printf("         compute the probability table!\n");
	      printf("  Building table for %d blocksize... ", nWordSize);
	      BuildChiTable( nWordSize, ctbl );
	      printf("\n");
	      if( nWordSize < 21) tblsz = nWordSize;
	      else tblsz = 21;
	      pvalue = longrun( &bitBlockltl, LITTLENDIAN, nWordSize, longRunBlocks, 
		ctbl, tblsz);
	      printf("Little endian p-value = %f\n", pvalue);
	      pvalue = longrun( &bitBlockbig, BIGENDIAN, nWordSize, longRunBlocks,
				ctbl, tblsz);
	      printf("Big endian p-value = %f\n", pvalue);
	    }
	    else
	    {
//		AnalyzeBitStream(fpInput, nWordSize, ii);
	    }
	  }
	  else
	  {
	    //	AnalyzeBitStream(fpInput, nWordSize, nBitNumber);
	  }
	}

	fclose(fpInput);
	
	return 0;
}
/*

  int FillBitStreamBuffer(BITSTREAM *bp)

  This function reads a pair of hex digits from the input data stream
  into a byte buffer.

  Function Returns:

  0 == Success.
  Anything else, some kind of I/O or data error

*/
int	FillBitStreamBuffer(BITSTREAM *bp)
{
	if (bp->nError == 0)
	{
		int		nFileChar;

		if (fscanf(bp->fp, " %2x", &nFileChar) == 1)
		{
			bp->ucByteBuf = (unsigned char) nFileChar;
		}
		else if ((bp->nError = ferror(bp->fp)) == 0)
		{
			bp->nError = (-1);
		}
	}

	return bp->nError;
}
/*

  int ResetBitStream(BITSTREAM *bp)

  This function resets the file pointer to the beginning
  of the file, reads the first byte into a buffer and

  Function returns:

  0 == Success
  Anything else, some kind of I/O error.

*/

int ResetBitStream(BITSTREAM *bp)
{
	if (bp->fp != NULL)
	{
		bp->nError = fseek(bp->fp, 0L, SEEK_SET);

		if (bp->nError == 0)
		{
			FillBitStreamBuffer(bp);
		}

		bp->nPrevBit = bp->nBitNumber - bp->nWordSize;
	}
	else
	{
		bp->nError = (-1);
	}

	return bp->nError;
}

/*

  int InitBitStream(BITSTREAM *bp, FILE *fp, int nWordSize)

  This function initializes a BITSTREAM object.

  Function returns:

  0 == Success
  Anything else, some kind of I/O error.

*/

int	InitBitStream(BITSTREAM *bp, FILE *fp, int nWordSize, int nBitNumber)
{
	bp->fp = fp;
	bp->nBitNumber = nBitNumber;
	bp->nWordSize = nWordSize;

	return ResetBitStream(bp);
}

/*

  int GetNextBit(BITSTREAM *bp)

  The function returns the next bit in a stream of bits
  from the specified input file.

  Function Returns:

  0  == The bit is zero
  1  == The bit is one
  -1 == Some kind of I/O or data error

*/

int GetNextBit(BITSTREAM *bp)
{
	int		nRet = (-1);

	if (bp->fp != NULL)
	{
		if (bp->nError == 0)
		{
			int		nBitNum = bp->nPrevBit + bp->nWordSize;
			int		nSkipBytes = nBitNum / BITS_PER_BYTE;
			int		nNextBit = nBitNum % BITS_PER_BYTE;

			while ((nSkipBytes > 0) && (bp->nError == 0))
			{
				if (FillBitStreamBuffer(bp) == 0)
				{
					nSkipBytes -= 1;
				}
			}

			if (bp->nError == 0)
			{
				nRet = (((int) bp->ucByteBuf & nByteMask[nNextBit]) == 0) ? 0 : 1;
				bp->nPrevBit = nNextBit;
			}
		}
	}

	return nRet;
}