⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 main.cpp

📁 随机数测试标准程序NIST
💻 CPP
📖 第 1 页 / 共 5 页
字号:
}

bool PokerTest(double &x)
{
	int f[16],i,j,k;
	double a=0,b=0;

	for(i=0; i<16; i++)
		f[i]=0;

	for(i=0; i<5000; i++)
	{
		j=i*4;
		k=8*bit[j] + 4*bit[j+1] + 2*bit[j+2] + bit[j+3];
		f[k]++;
	}

	for(i=0; i<16; i++)
	{
		a += pow(f[i], 2);
		//cout<<endl<<i <<"="<<f[i]<<endl;
		//b += f[i];
    }
    //cout<<endl<<"b="<<b<<endl;
    
	x = a*16/5000 - 5000;

	if(x>1.03 && x<57.4)
		return true;
	else
		return false;
}

bool RunsTest(int num[])
{
	int i,n=0;
	bool m;

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

	m=bit[0];
	for(i=0; i<20000; i++)
	{
		if(bit[i] == m)
			n++;
		else
		{
			if(n>6)
				n=6;
			num[n-1]++;
			m=bit[i];
			n=0;
		}
	}

	if(    num[0]>=2267 && num[0]<=2773 && num[1]>=1079 && num[1]<=1421
		&& num[2]>=502  && num[2]<=748  && num[3]>=223  && num[3]<=402
		&& num[4]>=90   && num[4]<=223  && num[5]>=90   && num[5]<=203)
		return true;
	else
		return false;
}

bool LongRunTest(int &max)
{
	int i,n=0;
	bool m;

	max=0;
	m=bit[0];
	for(i=0; i<20000; i++)
	{
		if(bit[i] == m)
			n++;
		else
		{
			if(n>max)
				max=n;
			m=bit[i];
			n=0;
		}
	}

	if(max < 34)
		return true;
	else
		return false;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                A P P R O X I M A T E  E N T R O P Y   T E S T
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

bool ApproximateEntropy(int m, int n)
{
	int           i, j, k, r, blockSize, seqLength;
	int           powLen, index;
	double        sum, numOfBlocks, ApEn[2], apen, chi_squared, p_value;
	unsigned int* P;
	

	seqLength = n;

	r = 0;

	for( blockSize=m; blockSize<=m+1; blockSize++ ) {
		if (blockSize == 0) {
			ApEn[0] = 0.00;
			r++;
		}
		else {
			numOfBlocks = (double)seqLength;
			powLen = (int)pow(2,blockSize+1)-1;
			if ( (P = (unsigned int*)calloc(powLen,sizeof(unsigned int)))== NULL ){
				printf("ApEn:  Insufficient memory available.\n");
				return false;
			}
			for( i=1; i<powLen-1; i++ )
				P[i] = 0;
			for( i=0; i<numOfBlocks; i++ ) { /* COMPUTE FREQUENCY */
				k = 1;
				for( j=0; j<blockSize; j++ ) {
					if ( (int)bit[(i+j)%seqLength] == 0 )
						k *= 2;
					else if ( (int)bit[(i+j)%seqLength] == 1 )
						k = 2*k+1;
				}
				P[k-1]++;
			}
			/* DISPLAY FREQUENCY */
			sum = 0.0;
			index = (int)pow(2,blockSize)-1;
			for( i=0; i<(int)pow(2,blockSize); i++ ) {
				if ( P[index] > 0 )
					sum += P[index]*log(P[index]/numOfBlocks);
				index++;
			}
			sum /= numOfBlocks;
			ApEn[r] = sum;
			r++;
			free(P);
		}
	}
	apen = ApEn[0] - ApEn[1];

	chi_squared = 2.0*seqLength*(log(2) - apen);
	p_value = (igamc)(pow(2,m-1),chi_squared/2.);

	if ( m > (int)(log(seqLength)/log(2)-5) ) {
		printf("\t\tNote: The blockSize = %d exceeds recommended", m);
		printf("\n\t\tvalues = %d; results are inaccurate!\n", MAX(1,(int)(log(seqLength)/log(2)-5)));
		printf("\t\t--------------------------------------------\n");
	}
	if ( p_value < ALPHA ) {
		return false;
	}
	else {
		return true;
	}
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                    B L O C K  F R E Q U E N C Y  T E S T
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

bool BlockFrequency(int m, int n)
{
	int    i, j, N;
	double blockSum, arg1, arg2, p_value;
	double sum, pi, v, chi_squared;
	
	N = (int)floor((double)n/(double)m); 		/* # OF SUBSTRING BLOCKS      */
	sum = 0.0;

	for( i=0; i<N; i++ ) {			/* N=10000 FOR EACH SUBSTRING BLOCK   */
		pi = 0.0;
		blockSum = 0.0;
		for( j=0; j<m; j++ ) 		/* m=100 COMPUTE The "i"th Pi Value */	 
			blockSum += bit[j+i*m];
		pi = (double)blockSum/(double)m;
		v = pi - 0.5;
		sum += v*v;
	}
	chi_squared = 4.0 * m * sum;

	arg1 = (double) N/2.e0;
	arg2 = chi_squared/2.e0;
	p_value = (igamc)(arg1,arg2);
	if ( p_value < ALPHA ) {
		return false;
	}
	else {
		return true;
	}
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
		    C U M U L A T I V E  S U M S  T E S T
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

bool CumulativeSums(int n)
{
	int    i, k, start, finish, mode;
	double p_value, cusum, z, sum, sum1, sum2;
		
	for( mode=0; mode<2; mode++ ) { /* mode = {0,1}  => {forward,reverse} */
		sum = 0.0;
		cusum = 1;

		if ( mode == 0 )
			for( i=0; i<n; i++ ) {
				sum += 2*(int)bit[i] - 1;
				cusum = MAX(cusum, fabs(sum));
			}
		else if (mode == 1)
			for( i=n-1; i>=0; i-- ) {
				sum += 2*(int)bit[i] - 1;
				cusum = MAX(cusum, fabs(sum));
			}
		z = cusum;

		sum1 = 0.0;
		start = (-n/(int)z+1)/4;
		finish = (n/(int)z-1)/4;
		for( k=start; k<=finish; k++ )
			sum1 += (normal((4*k+1)*z/sqrt(n))-normal((4*k-1)*z/sqrt(n)));

		sum2 = 0.0;
		start = (-n/(int)z-3)/4;
		finish = (n/(int)z-1)/4;
		for( k=start; k<=finish; k++ )
			sum2 += (normal((4*k+3)*z/sqrt(n))-normal((4*k+1)*z/sqrt(n)));
		p_value = 1.0 - sum1 + sum2;
		if ( p_value<0 || p_value>1)
			printf("\t\tWARNING:  P_VALUE IS OUT OF RANGE\n");
		if ( p_value < ALPHA ) {
			return false;
		}
		else {
			return true;
		}
	}
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
         D I S C R E T E  F O U R I E R  T R A N S F O R M  T E S T 
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
//void  __ogg_fdrffti(int n, double* wsave, double* ifac);
//void  __ogg_fdrfftf(int n, double* X, double* wsave, double* ifac);
bool DiscreteFourierTransform(int n)
{
	double   p_value, upperBound, percentile, *m, *X;
	int      i, count;
	double   N_l, N_o, d;

	double*  wsave, *ifac;
	  
	if ( ((X = (double*) calloc(n,sizeof(double))) == NULL) ||
			 ((wsave = (double *)calloc(2*n+15,sizeof(double))) == NULL) ||
			 ((ifac = (double *)calloc(15,sizeof(double))) == NULL) ||
			 ((m = (double*)calloc(n/2+1, sizeof(double))) == NULL) ) {
			printf("\t\tUnable to allocate working arrays for the DFT.\n");
			if( X != NULL )
				free(X);
			if( wsave != NULL )
				free(wsave);
			if( ifac != NULL )
				free(ifac);
			if( m != NULL )
				free(m);
		}
		else {
			for( i=0; i<n; i++ )
				X[i] = 2*(int)bit[i] - 1;

			__ogg_fdrffti(n, wsave, ifac);		    /* INITIALIZE WORK ARRAYS */
			__ogg_fdrfftf(n, X, wsave, ifac);	    /* APPLY FORWARD FFT */

			m[0] = sqrt(X[0]*X[0]);	    /* COMPUTE MAGNITUDE */

			for( i=0; i<n/2; i++ ) {	   	    /* DISPLAY FOURIER POINTS */
				m[i+1] = sqrt( pow(X[2*i+1],2) + pow(X[2*i+2],2) ); 
			}
			count = 0;				       /* CONFIDENCE INTERVAL */
			//upperBound = sqrt(3*n);
			upperBound = sqrt(2.995732274*n);
			for( i=0; i<n/2; i++ )
				if ( m[i] < upperBound )
					count++;
			percentile = (double)count/(n/2)*100;
			N_l = (double) count;       /* number of peaks less than h = sqrt(3*n) */
			N_o = (double) 0.95*n/2.;
			//d = (N_l - N_o)/sqrt(n/2.*0.95*0.05);
			d = (N_l - N_o)/sqrt(n/4.0*0.95*0.05);
			p_value = erfc(fabs(d)/sqrt(2.));
			free(X);
			free(wsave);
			free(ifac);
			free(m);
			if ( p_value < ALPHA ) {
				return false;
			}
			else {
				return true;
			}
		}
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                          F R E Q U E N C Y  T E S T
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

bool Frequency(int n)
{
	int    i;
	double f, s_obs, p_value, sum;
	double sqrt2 = 1.41421356237309504880;
		
	sum = 0.0;
	for( i=0; i<n; i++ )
		sum += 2*(int)bit[i]-1;
	s_obs = fabs(sum)/sqrt(n);
	f = s_obs/sqrt2;
	p_value = erfc(f);
	if ( p_value < ALPHA ) {
			return false;
	}
	else {
			return true;
	}
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
               L E M P E L  Z I V  C O M P R E S S I O N  T E S T
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

bool LempelZivCompression(int n)
{
	int			W;	/* Number of words */
	int			i, j, k, prev_I, powLen, max_len, nOrig;
	int			done = 0;
	double		p_value, mean, variance; 
	BitField	*P;
	
	W = 0;
	nOrig = n;
	k = (int)(log(n)/log(2)+6);
	powLen = (int)pow(2, k);
	if ( (P = (BitField*)calloc(powLen, sizeof(BitField))) == NULL ) {
		printf("\t\tLempel-Ziv Compression has been aborted,\n");    
		printf("\t\tdue to insufficient workspace!\n");
	}
	else {
		for( i=0; i<powLen; i++ ) 
			P[i].b = 0;
		i = 0;
		max_len = 0;
		while( i <= n-1 ) {
			done = 0;
			j = 1;
			prev_I = i;
			while( !done ) {
				if ( 2*j+1 <= powLen ) {
					if ( (int)bit[i] == 0 ) {
						if ( P[2*j].b == 1 ) {
							j *= 2;
						}
						else {
							P[2*j].b = 1;
							done = 1;
						}
					}
					else {
						if ( P[2*j+1].b == 1 ) {
							j = 2*j+1;
						}
						else {
							P[2*j+1].b = 1;
							done = 1;
						}
					}
					i++;
					if ( i > n-1 ) {
						done = 1;
					}
					if ( done ) {
						W++;
						max_len = MAX(max_len, i-prev_I);
					}
				}
				else {
					printf("\t\tWARNING: Segmentation Violation Imminent.");
					printf("\n\t\t\t Lempel-Ziv Compression Terminated.\n");
					printf("\t\t-----------------------------------------\n");
					done = 1;
					i = n;
				}
			}
		}
		switch( n ) {
		case 100000:
			mean = 8782.500000; 
			variance = 11.589744;
			break;
		case 200000:
			mean = 16292.1000;
			variance = 21.4632;
			break;
		case 400000:
			mean = 30361.9500;
			variance = 58.7868;
			break;
		case 600000:
			mean = 43787.5000;
			variance = 70.1579;
			break;
		case 800000:
			mean = 56821.9500;
			variance = 67.4184;
			break;
		case 1000000:
			/* Updated Mean and Variance 10/26/99 */
			mean = 69588.20190000;
			variance = 73.23726011;
			/* Previous Mean and Variance
			mean = 69586.250000;
			variance = 70.448718;
			*/
			break;
		case 5000000:
			mean = 303156.9;
			variance = 266.1341;
			//fp = fopen("W-obs-5000000.txt", "a");
			//fprintf(fp, "Wobs: %d\n", W);
			//fclose(fp);
			break;
		case 10000000:
			mean = 574227.482;
			variance = 447.2809059;
			//fp = fopen("W-obs.txt", "a");
			//fprintf(fp, "Wobs: %d\n", W);
			//fclose(fp);
			break;
		case 100000000:
			mean = 574227.482;
			variance = 447.2809059;
//			fp1 = fopen("W-obs-100000000.txt", "a");
//			fprintf(fp1, "%d\n", W);
//			fclose(fp1);
			break;
		default:
			break;
		}
		p_value = 0.5*erfc((mean-W)/pow(2.0*variance,0.5));
		free(P);
		if ( p_value < ALPHA ) {
				return false;
		}

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -