main.cpp

随机数测试标准程序NIST

#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <float.h>
#include "proto.h"
#include "special-functions.h"
#include "cephes-protos.h"
#define  ALPHA 0.01
#include "defs.h"
#include "externs.h"
#include "dfft.h"
#include "utilities1.h"
#include "matrix.h"
using namespace std;

bool bit[20000],buf[32];

static int p = 0, point = 0;

bool Von_adjust(char dat, bool *bit);
bool SHA_1(bool *buf, bool *bit);
void       perform_elementary_row_operations(int,int,int,int,BitField**);
int        find_unit_element_and_swap(int,int,int,int, BitField**);
int        swap_rows(int, int, int, BitField**);
int        determine_rank(int, int, int, BitField**);
int        computeRank(int,int,BitField**);
void       define_matrix(int,int,int,BitField**);
BitField** create_matrix(int, int);
void       display_matrix(int, int, BitField**);
void       def_matrix(int, int, BitField**, int);
void       delete_matrix(int, BitField**);

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                 R A N K  A L G O R I T H M  R O U T I N E S
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

int computeRank(int M, int Q, BitField** matrix)
{
   int i;
   int rank;
   int m = MIN(M,Q);

   /* FORWARD APPLICATION OF ELEMENTARY ROW OPERATIONS */ 
      
   for(i = 0; i < m-1; i++) {
      if (matrix[i][i].b == 1) 
         perform_elementary_row_operations(0, i, M, Q, matrix);
      else { 	/* matrix[i][i] = 0 */
         if (find_unit_element_and_swap(0, i, M, Q, matrix) == 1) 
            perform_elementary_row_operations(0, i, M, Q, matrix);
      }
      if (MATRICES) {
         //fprintf(output,"FORWARD ITERATION = %d\n", i);
         display_matrix(M,Q, matrix);
      }
   }

   /* BACKWARD APPLICATION OF ELEMENTARY ROW OPERATIONS */ 
   for(i = m-1; i > 0; i--) {
      if (matrix[i][i].b == 1)
         perform_elementary_row_operations(1, i, M, Q, matrix);
      else { 	/* matrix[i][i] = 0 */
         if (find_unit_element_and_swap(1, i, M, Q, matrix) == 1) 
            perform_elementary_row_operations(1, i, M, Q, matrix);
      }
      if (MATRICES) {
        // fprintf(output,"BACKWARD ITERATION = %d\n", i);
         display_matrix(M,Q, matrix);
      }
   } 
   #ifdef BACK
  // fprintf(output, "\n");
   display_matrix(M,Q, matrix);
   #endif
   rank = determine_rank(m,M,Q,matrix);
  // if (MATRICES) fprintf(output,"RANK = %d\n", rank);
   return rank;
}

void perform_elementary_row_operations(int flag, int i, int M, int Q, 
			               BitField** A)
{
  int j,k;

  switch(flag)
  { 
     case 0: for(j = i+1; j < M;  j++)
                if (A[j][i].b == 1) 
                   for(k = i; k < Q; k++) 
                      A[j][k].b = (A[j][k].b + A[i][k].b) % 2;
             break;
        
     case 1: for(j = i-1; j >= 0;  j--)
                if (A[j][i].b == 1)
                   for(k = 0; k < Q; k++)
                       A[j][k].b = (A[j][k].b + A[i][k].b) % 2;
             break;
     default:  printf("ERROR IN CALL TO perform_elementary_row_operations\n");
               break;
  }
  return;
}

int find_unit_element_and_swap(int flag, int i, int M, int Q, BitField** A)
{ 
  int index;
  int row_op = 0;

  switch(flag) 
  {
    case 0:  index = i+1;
             while ((index < M) && (A[index][i].b == 0)) 
                index++;
             if (index < M) 
                row_op = swap_rows(i,index,Q,A);
             break;
    case 1:
             index = i-1;
	     while ((index >= 0) && (A[index][i].b == 0)) 
	       index--;
	     if (index >= 0) 
                row_op = swap_rows(i,index,Q,A);
             break;
    default:  printf("ERROW IN CALL TO find_unit_element_and_swap\n");
              break;
  }
  return row_op;
}

int swap_rows(int i, int index, int Q, BitField** A)
{
  int p;
  BitField temp;

  for(p = 0; p < Q; p++) {
     temp.b = A[i][p].b;
     A[i][p].b = A[index][p].b;
     A[index][p].b = temp.b;
  }
  return 1;
}

int determine_rank(int m, int M, int Q, BitField** A)
{
   int i, j, rank, allZeroes;

   /* DETERMINE RANK, THAT IS, COUNT THE NUMBER OF NONZERO ROWS */

   rank = m;
   for(i = 0; i < M; i++) {
      allZeroes = 1; 
      for(j=0; j < Q; j++) {
         if (A[i][j].b == 1) {
            allZeroes = 0;
            break;
         }
      }
      if (allZeroes == 1) rank--;
   } 
   return rank;
}

BitField** create_matrix(int M, int Q)
{
  int i;
  BitField **matrix;

  if ((matrix = (BitField**) calloc(M, sizeof(BitField *))) == NULL) {
     printf("ERROR IN FUNCTION create_matrix:  Insufficient memory available.");
     printf("\n");
     printf("CREATE_MATRIX:  Insufficient memory for %dx%d matrix.\n",
	     M, M); 
     return matrix;	
  }
  else {
     for (i = 0; i < M; i++) {
        if ((matrix[i] = (BitField*)calloc(Q, sizeof(BitField))) == NULL) {
           printf("CREATE_MATRIX:  Insufficient memory for %dx%d ",M,M);
           printf("matrix.\n"); 
           printf("ERROR IN FUNCTION create_matrix: Insufficient memory for ");
           printf("%dx%d matrix.\n", M, M); 
	   return NULL;
        }
     }
     return matrix;
  }
}

void display_matrix(int M, int Q, BitField** m)
{
  int i,j;

  for (i = 0; i < M; i++) {
     for (j = 0; j < Q; j++)
        printf("%d ", m[i][j].b);
     printf("\n");
  }
  return;
}

void def_matrix(int M, int Q, BitField** m,int k)
{
  int   i,j;
  for (i = 0; i < M; i++) 
     for (j = 0; j < Q; j++) { 
         m[i][j].b = bit[k*(M*Q)+j+i*M];
     }
  return;
}

void delete_matrix(int M, BitField** matrix)
{
  int i;
  for (i = 0; i < M; i++)
    free(matrix[i]);
  free(matrix);
}


//extern double erf ( double x );

double normal(double x)
{
	double arg, result, sqrt2=1.414213562373095048801688724209698078569672;

	if (x > 0) {
		arg = x/sqrt2;
		result = 0.5 * ( 1 + erf(arg) );
	}
	else {
		arg = -x/sqrt2;
		result = 0.5 * ( 1 - erf(arg) );
	}

	return(result);
}

double normal2(double a)
{
	return (1.0-0.5*erfc(a/sqrt(2.0)));
}

int Getdata(void)
{
    char ch;
    unsigned int dat;
    ifstream data("data.txt",ios::in);
    if(!data)
    {
      cout<<"open data.txt fail!"<<endl;
      return (-1);
    }
    else
    {
      data.seekg(0,ios::beg);
      while(!data.eof() && (point<20000))
      {
        data.get(ch);
        if((ch != '\n')&&(ch != '\r') )
        {
             if(ch == ' ')
             {
                ch = '0';
             }
            dat = atoi(&ch);
            if((ch>='a')&&(ch<='f'))
             dat = ch - 'a' + 10;
            Von_adjust( dat, bit);
        }
       // else
        //cout<<endl;//hex<<dat<<endl;//<<"ch="<<hex<<ch<<endl;
      }
      if(point<20000)
        cout<<"data < 20000!!!"<<endl;
      //else
        //cout<<"Get data!"<<endl;
      data.clear();
    }
    return (0);
}

bool Von_adjust(char dat, bool *bit)
{
     if((dat>>3) != (dat&4)>>2)
     {
         if((dat>>3) == 1)
             buf[p++] = 0;
         else
             buf[p++] = 1;
         if( p == 32)
         {
             p = 0;
             if(SHA_1(buf, bit))
               return 1;
         }         
     }
     if((dat&2)>>1 != (dat&1))
     {
         if((dat&2)>>1 == 1 )
             buf[p++] = 0;
         else
             buf[p++] = 1; 
         if( p == 32)
         {
             p = 0;
             if(SHA_1(buf, bit))
               return 1;
         }
     }
     return 0;
}

unsigned int f(unsigned int t,unsigned int x, unsigned int y, unsigned int z)
{
    if(t>=0 && t<=19 )
        return ((x&y)|((~x)&z));
    else if(t>=20 && t<=39)
        return (x^y^z);
    else if(t>=40 && t<=59)
        return ((x&y)|(x&z)|(y&z));
    else if(t>=60 && t<=79)
        return (x^y^z);
	else
		return 0;
}

unsigned int K(unsigned int t)
{
    if(t>=0 && t<=19)
         return 0x5a827999;
    else if(t>=20 && t<=39)
         return 0x6ed9eba1;
    else if(t>=40 && t<=59)
         return 0x8f1bbcdc;
    else if(t>=60 && t<=79)
         return 0xca62c1d6;
	else
		 return 0;
}

bool SHA_1(bool *buf, bool *bit)
{
     //bool reg[512];
     unsigned int w[80],a,b,c,d,e,T,H0,H1,H2,H3,H4;
     int i;
     
     H0 = 0x67452301;
     H1 = 0xefcdab89;
     H2 = 0x98badcfe;
     H3 = 0x10325476;
     H4 = 0xc3d2e1f0;

     for(i=0; i<80; i++)
     {
         if( i<4 )
             for(int j=0; j<8; j++)
                 w[i] |= buf[8*i+j]<<(7-j);
         else if( i==5 )
             w[i] = 0x80;
         else if( i<15 )
             w[i] = 0;
         else if( i==15 )
             w[i] = 24;//32;
         else if( i>15 )
             w[i] = ((w[i-3]^w[i-8]^w[i-14]^w[i-16])<<1)|((w[i-3]^w[i-8]^w[i-14]^w[i-16])>>31);
     }
     a = H0;
     b = H1;
     c = H2;
     d = H3;
     e = H4; 
     for(i=0; i<80; i++)
     {
        T = ((a<<5)|(a>>27))+f(i,b,c,d)+e+K(i)+w[i];
        e = d;
        d = c;
        c = ((b<<30)|(b>>2));
        b = a;
        a = T;
#if 0
        if(i<20||i>76)
          {
            cout<<dec<<i<<endl;
            cout<<hex<<K(i)<<endl;
            cout<<hex<<a<<endl<<b<<endl<<c<<endl<<d<<endl<<e<<endl<<endl;
          }
#endif          
     }
     H0 = a+H0;
     H1 = b+H1;
     H2 = c+H2;
     H3 = d+H3;
     H4 = e+H4;
     //cout<<hex<<H0<<endl<<H1<<endl<<H2<<endl<<H3<<endl<<H4<<endl;
     for(i=0; i<32; i++)
     {
        if(((H0<<i)>>31) == 1)
            bit[point++] = 1;
        else
            bit[point++] = 0;
        if(point == 20000)
            return 0; 
     }
     for(i=0; i<32; i++)
     {
        if(((H1<<i)>>31) == 1)
            bit[point++] = 1;
        else
            bit[point++] = 0;
        if(point == 20000)
            return 0; 
     }
     for(i=0; i<32; i++)
     {
        if(((H2<<i)>>31) == 1)
            bit[point++] = 1;
        else
            bit[point++] = 0;
        if(point == 20000)
            return 0; 
     }
     for(i=0; i<32; i++)
     {
        if(((H3<<i)>>31) == 1)
            bit[point++] = 1;
        else
            bit[point++] = 0;
        if(point == 20000)
            return 0; 
     }
     for(i=0; i<32; i++)
     {
        if(((H4<<i)>>31) == 1)
            bit[point++] = 1;
        else
            bit[point++] = 0;
        if(point == 20000)
            return 0; 
     }
     return 1;
}

bool MonobitTest(int &num)
{
	num=0;
	for(int i=0; i<20000; i++)
		if(bit[i])
			num++;

	if(num>9654 && num <10364)
		return true;
	else
		return false;