МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ ,,ЛЬВІВСЬКА ПОЛІТЕХНІКА’’
Кафедра ІСМ
Лабораторна робота № 2
з дисципліни “Моделювання систем”
на тему
“Оцінювання якості генераторів псевдовипадкових чисел”
Львів-2013
Мета роботи
Вивчення критеріїв якості генераторів випадкових чисел та програмна реалізація оцінювання якості послідовностей псевдовипадкових чисел.
Варіант завдання
Варіант 17
Код основної програми
#include <iostream>
#include <math.h>
#include <windows.h>
#include <stdint.h>
#include <ctime>
#include <cstdlib>
#include <string>
#include <sstream>
using namespace std;
#define m 1000000
#pragma comment(linker, "/STACK:100000000")
#define PI 3.14159265358979323846
double result;
double XSK=0.37470302; //Х0 середніх квадратів
int g=16;
double XMM=28768; //Х0 мультиплікативного методу
int t=162;
int l=811;
int Seed=8005;
unsigned Array[624]; //массив ПВЧ
unsigned Y; //змінна, що буде містити ПВЧ у вихорі Мерсенна
int index; //індекс, який буде "бігати" по масиву
double b;
//ТАЙМЕР
class CTimer
{
public:
CTimer()
{
QueryPerformanceFrequency(&mqFreq);
}
~CTimer() {}
void Start()
{
QueryPerformanceCounter(&mqStart);
}
void End()
{
QueryPerformanceCounter(&mqEnd);
}
double GetTimeInSeconds()
{
return (mqEnd.QuadPart - mqStart.QuadPart)/static_cast<double>(mqFreq.QuadPart);
}
private:
LARGE_INTEGER mqStart;
LARGE_INTEGER mqEnd;
LARGE_INTEGER mqFreq;
};
//МЕТОД СЕРЕДИННИХ КВАДРАТІВ
double MSK()
{
double temp; //тимчасова змінна
double *P=0; //вказівник на ділянку пам'яті, в якій буде зберігатись ціла частина числа
P=new double(sizeof(double)); //виділення пам'яті під вказівник
double XSKtemp=XSK*XSK; //Х^2
XSKtemp=XSKtemp*1000000000000; //|
temp=modf(XSKtemp,P); //|позбавляємось лівих бітів
XSKtemp=*P; //|
XSKtemp=XSKtemp/100000000; //|
temp=modf(XSKtemp,P); //|позбавляємось правих бітів
XSKtemp=temp;
XSK=XSKtemp; //змінюю значення Х0
delete P;
return XSKtemp;
}
//МУЛЬТИПЛІКАТИВНИЙ МЕТОД
double Multiplicative_Method()
{
uint16_t X1=XMM*l; //беремо 16 молодших бітів числа
XMM=X1; //змінюю значення Х0
result=X1;
result/=65536; //переводжу випадкове чесле до інтервалу (0,1)
return result;
}
double Mersenns_Whirl()
{
Y=Array[index]; //|
Y=Y^(Y>>11); //|
Y=Y^((Y<<7)&(2636928640)); //| генерація ПВЧ
Y=Y^((Y<<15)&(4022730752)); //|
Y=Y^(Y>>18); //|
index++;
result=Y;
result/=4294967296;
if(index>=624)
{
for(int j=0;j<624;j++) //генерація нових чисел для масиву
{
unsigned Z=(Array[j]&(0x80000000))+((Array[(j+1)%624])&(0x7FFFFFFF));
if((Z%2)==0)
{
Array[j]=Array[(j+397)%624]^(Z>>1);
}
else
{
Array[j]=Array[(j+397)%624]^(Z>>1)^(2567483615);
}
}
index=0; //повертаємось до початку масиву
}
return result;
}
double Random()
{
double RandDigit;
RandDigit=rand()/(double)RAND_MAX;
return RandDigit;
}
double Gistogram(double Digits_Array[], int intervals)
{
double deviation=0;
double Array1[10];
double Array2[100];
double Array3[1000];
double PVCh10;
if(intervals==10)
{
for(int i=0; i<10; i++)
Array1[i]=0;
for(int i=0; i<m; i++)
{
PVCh10=Digits_Array[i]*10;
int flag=PVCh10; //flag містить номер інтервалу, у який попало число
*(Array1+flag)=*(Array1+flag)+1;
}
/*for(int i=0; i<10; i++)
cout<<Array1[i]<<endl; */ //кількість чисел в кожному інтервалі
for(int i=0; i<intervals; i++)
Array1[i]-=100000;
for(int i=0; i<intervals; i++)
Array1[i]*=Array1[i];
for(int i=0; i<intervals; i++)
deviation+=Array1[i];
deviation/=intervals;
}
/////////////////////////////////////////////////////////////
if(intervals==100)
{
for(int i=0; i<100; i++)
Array2[i]=0;
for(int i=0; i<m; i++)
{
PVCh10=Digits_Array[i]*100;
int flag=PVCh10; //flag містить номер інтервалу, у який попало число
*(Array2+flag)=*(Array2+flag)+1;
}
/*for(int i=0; i<100; i++)
cout<<Array2[i]<<endl;*/
for(int i=0; i<intervals; i++)
Array2[i]-=10000;
for(int i=0; i<intervals; i++)
Array2[i]*=Array2[i];
for(int i=0; i<intervals; i++)
deviation+=Array2[i];
deviation/=intervals;
}
//////////////////////////////////////////////////
if(intervals==1000)
{
for(int i=0; i<1000; i++)
Array3[i]=0;
for(int i=0; i<m; i++)
{
PVCh10=Digits_Array[i]*1000;
int flag=PVCh10; //flag містить номер інтервалу, у який попало число
*(Array3+flag)=*(Array3+flag)+1;
}
/*for(int i=0; i<100; i++)
cout<<Array3[i]<<endl;*/
for(int i=0; i<intervals; i++)
Array3[i]-=1000;
for(int i=0; i<intervals; i++)
Array3[i]*=Array3[i];
for(int i=0; i<intervals; i++)
deviation+=Array3[i];
deviation/=intervals;
}
return deviation;
}
double Ring(double Digits_Array[])
{
int ring=0;
for(int i=0;i<m;i+=2)
{
if((pow(Digits_Array[i],2)+pow(Digits_Array[i+1],2))<1)
ring++;
}
cout<<ring;
double ideal=(m/2.0)*(PI/4.0);
result=ring/ideal;
return result;
}
long double factorial(long int d)//факторіал по іншому рахувати треба!
{
long double res=1;
for (int i=2;i<=d;i++)
res*=i;
/*cout<<res<<endl;
system("pause");*/
return res;
}
double Stochasticity(double Digits_Array[])
{
int l=32;
double temp;
double Pjl[33];
int counter[33];
int ideal_counter[33];
for(int i=0;i<33;i++)
counter[i]=0;
unsigned int digit;
int shift=0;
for(int i=0;i<m;i++)
{
digit=Digits_Array[i]*4294967296;
digit=digit-((digit>>1)&0x55555555);
digit=(digit& 0x33333333)+((digit>>2)&0x33333333);
digit=(digit+(digit>>4))&0x0F0F0F0F;
digit=digit+(digit>>8);
digit=digit+(digit>>16);
shift=digit&0x0000003F;
*(counter+shift)=*(counter+shift)+1;
}
/*for(int i=0;i<33;i++)
cout<<i<<":"<<counter[i]<<endl;*/
long double cjl;
for(int i=0;i<33;i++)
{
cjl=factorial(l)/(factorial(i)*factorial(l-i));
temp=pow(0.5,l);
Pjl[i]=temp*cjl;
}
/*
double PP=0;
for(int i=0;i<33;i++)
PP+=Pjl[i];
cout<<"Сума iмовiрностей = "<<PP<<endl;
*/
for(int i=0;i<33;i++)
{
ideal_counter[i]=Pjl[i]*1000000;
}
result=0;
for(int i=0;i<33;i++)
result+=pow((double)(counter[i]-ideal_counter[i]),2);
result/=33;
return result;
}
int Period_length(double Digits_Array[])
{
int length;
for(int i=m-2;i>=0;i--)
if(Digits_Array[m-1]==Digits_Array[i])
{length=m-i; break;}
else
length=0;
return length;
}
double Independence(double Digits_Array[], int t)
{
double k1=0,k2=0,k3=0;
double Dxi,Dxit;
/////////////////////////
double tau=(double)t;
for (int i=0;i<m-t-1;i++)
k1+=Digits_Array[i];
k1=k1*(1/(m-tau));
for (int j=0;j<m-t-1;j++)
k2+=Digits_Array[j];
k2*=k2;
k2*=(1/((m-tau)*(m-tau)));
k1-=k2;
Dxi=k1;
/////////////////////////
k1=0; k2=0;
for (int i=0;i<m-t-1;i++)
k1+=Digits_Array[i+t];
k1=k1*(1/(m-tau));
for (int j=0;j<m-t-1;j++)
k2+=Digits_Array[j+t];
k2*=k2;
k2*=(1/((m-tau)*(m-tau)));
k1-=k2;
Dxit=k1;
/////////////////////////
k1=0; k2=0;
for (int i=0;i<m-t-1;i++)
k1+=Digits_Array[i]*Digits_Array[i+t];
k1=k1*(1/(m-tau));
for (int j=0;j<m-t-1;j++)
k2+=Digits_Array[j+t];
for (int l=0;l<m-t-1;l++)
k3+=Digits_Array[l];
k2*=k3;
k2=k2*(1/((m-tau)*(m-tau)));
k1-=k2;
k3=Dxi*Dxit;
result=k1/pow((double)k3,0.5);
b=result/pow((1/m),0.5);
return result;
}
//ГОЛОВНА ФУНКЦІЯ
void main()
{
double Digits_Array[m];
int intervals=0;
Array[0]=Seed; //початкове значення
for(int j=1;j<624;j++)
{
Array[j]=1812433253*(Array[j-1]^(Array[j-1]>>30))+j; //ініціалізація масиву
}
setlocale(LC_ALL, "Ukrainian");
CTimer Time;
cout<<"Кiлькiсть iнтервалiв(10/100/1000): ";
for(;;)
{
cin>>intervals;
if((intervals==10) || (intervals==100) || (intervals==1000))
break;
}
//Заповнення масиву випадковими числами за методом середніх квадратів
cout<<endl<<"Метод середнiх квадратiв"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=MSK();
cout<<"Середньо квадратична похибка = "<<Gistogram(Digits_Array, intervals)<<endl;
cout<<endl;
//Заповнення масиву випадковими числами за мультиплікативним методом
cout<<"Мультиплiкативни метод"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Multiplicative_Method();
cout<<"Середньо квадратична похибка = "<<Gistogram(Digits_Array, intervals)<<endl;
cout<<endl;
//Заповнення масиву випадковими числами за методом вихора Мерсенна
cout<<"Метод вихора Мерсенна"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Mersenns_Whirl();
cout<<"Середньо квадратична похибка = "<<Gistogram(Digits_Array, intervals)<<endl;
cout<<endl;
//Заповнення масиву випадковими числами з використанням стандартного генератора
cout<<"Стандартний генератор"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Random();
cout<<"Середньо квадратична похибка = "<<Gistogram(Digits_Array, intervals)<<endl;
cout<<endl<<"*******************************"<<endl;
XSK=0.21878104;
cout<<"Перевiрка за непрямими ознаками"<<endl<<"*******************************"<<endl;
cout<<"Метод середнiх квадратiв"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=MSK();
cout<<"Вiдношення кiлькоcтi попадань до теоретично очiкованої = "<<Ring(Digits_Array)<<endl<<endl;
XMM=30121;
t=101;
l=811;
cout<<"Мультиплiкативний метод"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Multiplicative_Method();
cout<<"Вiдношення кiлькоcтi попадань до теоретично очiкованої = "<<Ring(Digits_Array)<<endl<<endl;
Seed=1193;
Array[0]=Seed; //початкове значення
for(int j=1;j<624;j++)
{
Array[j]=1812433253*(Array[j-1]^(Array[j-1]>>30))+j; //ініціалізація масиву
}
cout<<"Метод вихора Мерсенна"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Mersenns_Whirl();
cout<<"Вiдношення кiлькоcтi попадань до теоретично очiкованої = "<<Ring(Digits_Array)<<endl<<endl;
cout<<"Стандартний генератор"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Random();
cout<<"Вiдношення кiлькоcтi попадань до теоретично очiкованої = "<<Ring(Digits_Array)<<endl<<endl;
cout<<endl<<"*******************************"<<endl;
cout<<"Стохастичнiсть"<<endl<<"*******************************"<<endl;
XSK=0.21878104;
cout<<endl<<"Метод середнiх квадратiв"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=MSK();
cout<<"Середньоквадратична похибка ймовiрностей появи рiзних кiлькостей бiтiв-одиниць "<<endl<<Stochasticity(Digits_Array)<<endl;
XMM=30121;
t=101;
l=811;
cout<<endl<<"Мультиплiкативний метод"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Multiplicative_Method();
cout<<"Середньоквадратична похибка ймовiрностей появи рiзних кiлькостей бiтiв-одиниць "<<endl<<Stochasticity(Digits_Array)<<endl;
Seed=1193;
Array[0]=Seed; //початкове значення
for(int j=1;j<624;j++)
{
Array[j]=1812433253*(Array[j-1]^(Array[j-1]>>30))+j; //ініціалізація масиву
}
cout<<endl<<"Метод вихора Мерсенна"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Mersenns_Whirl();
cout<<"Середньоквадратична похибка ймовiрностей появи рiзних кiлькостей бiтiв-одиниць "<<endl<<Stochasticity(Digits_Array)<<endl;
cout<<endl<<"Стандартний генератор"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Random();
cout<<"Середньоквадратична похибка ймовiрностей появи рiзних кiлькостей бiтiв-одиниць "<<endl<<Stochasticity(Digits_Array)<<endl;
XSK=0.21878104;
cout<<endl;
cout<<endl<<"Метод середнiх квадратiв"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=MSK();
cout<<"Довжина перiоду = "<<Period_length(Digits_Array)<<endl;
XMM=30121;
t=101;
l=811;
cout<<endl<<"Мультиплiкативний метод"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Multiplicative_Method();
cout<<"Довжина перiоду = "<<Period_length(Digits_Array)<<endl;
Seed=1193;
Array[0]=Seed; //початкове значення
for(int j=1;j<624;j++)
{
Array[j]=1812433253*(Array[j-1]^(Array[j-1]>>30))+j; //ініціалізація масиву
}
cout<<endl<<"Метод вихора Мерсенна"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Mersenns_Whirl();
cout<<"Довжина перiоду = "<<Period_length(Digits_Array)<<endl;
cout<<endl<<"Стандартний генератор"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Random();
cout<<"Довжина перiоду = "<<Period_length(Digits_Array)<<endl;
int t1=1;
XSK=0.21878104;
cout<<endl;
cout<<endl<<"Метод середнiх квадратiв"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=MSK();
double t=Independence(Digits_Array,1);
for(int i=2;i<=100;i++)
if(Independence(Digits_Array,i)>t){t1=i;}
cout<<"Кореляцiйний момент = "<<Independence(Digits_Array,t)<<" t = "<<t1<<endl;
t=0; t1=1;
XMM=30121;
t=101;
l=811;
cout<<endl<<"Мультиплiкативний метод"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Multiplicative_Method();
t=Independence(Digits_Array,1);
for(int i=2;i<=100;i++)
if(Independence(Digits_Array,i)>t){t1=i;}
cout<<"Кореляцiйний момент = "<<Independence(Digits_Array,100)<<" t = "<<t1<<endl;
t=0; t1=1;
Seed=1193;
Array[0]=Seed; //початкове значення
for(int j=1;j<624;j++)
{
Array[j]=1812433253*(Array[j-1]^(Array[j-1]>>30))+j; //ініціалізація масиву
}
cout<<endl<<"Метод вихора Мерсенна"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Mersenns_Whirl();
t=Independence(Digits_Array,1);
for(int i=2;i<=100;i++)
if(Independence(Digits_Array,i)>t){t1=i;}
cout<<"Кореляцiйний момент = "<<Independence(Digits_Array,100)<<" t = "<<t1<<endl;
t=0; t1=1;
cout<<endl<<"Стандартний генератор"<<endl;
for(int i=0;i<m;i++)
Digits_Array[i]=Random();
t=Independence(Digits_Array,1);
for(int i=2;i<=100;i++)
if(Independence(Digits_Array,i)>t){t1=i;}
cout<<"Кореляцiйний момент = "<<Independence(Digits_Array,100)<<" t = "<<t1<<endl;
system("pause");
}
Результат роботи програми:
Кiлькiсть iнтервалiв(10/100/1000): 10
Метод середнiх квадратiв
Середньо квадратична похибка = 8.9902e+010
Мультиплiкативни метод
Середньо квадратична похибка = 2.13487e+006
Метод вихора Мерсенна
Середньо квадратична похибка = 62552.4
Стандартний генератор
Середньо квадратична похибка = 54618.7
*******************************
Перевiрка за непрямими ознаками
*******************************
Метод середнiх квадратiв
387711Вiдношення кiлькоcтi попадань до теоретично очiкованої = 0.987799
Мультиплiкативний метод
392994Вiдношення кiлькоcтi попадань до теоретично очiкованої = 1.00126
Метод вихора Мерсенна
392739Вiдношення кiлькоcтi попадань до теоретично очiкованої = 1.00061
Стандартний генератор
392127Вiдношення кiлькоcтi попадань до теоретично очiкованої = 0.99905
*******************************
Стохастичнiсть
*******************************
Метод середнiх квадратiв
Середньоквадратична похибка ймовiрностей появи рiзних кiлькостей бiтiв-одиниць
5.27579e+006
Мультиплiкативний метод
Середньоквадратична похибка ймовiрностей появи рiзних кiлькостей бiтiв-одиниць
7.09849e+009
Метод вихора Мерсенна
Середньоквадратична похибка ймовiрностей появи рiзних кiлькостей бiтiв-одиниць
29047.9
Стандартний генератор
Середньоквадратична похибка ймовiрностей появи рiзних кiлькостей бiтiв-одиниць
3.32218e+008
Метод середнiх квадратiв
Довжина перiоду = 3166
Мультиплiкативний метод
Довжина перiоду = 16385
Метод вихора Мерсенна
Довжина перiоду = 0
Стандартний генератор
Довжина перiоду = 5075
Метод середнiх квадратiв
Кореляцiйний момент = 0.331965 t = 80
Мультиплiкативний метод
Кореляцiйний момент = -0.000466021 t = 91
Метод вихора Мерсенна
Кореляцiйний момент = 0.000221623 t = 100
Стандартний генератор
Кореляцiйний момент = 0.000359566 t = 70
Для продолжения нажмите любую клавишу . . .
Гістограми
Метод середніх квадратів
Мультиплікативний метод
Метод вихора Мерсенна
Вбудований генератор
Середньоквадратичні похибки
Оцінки попадань в одиничне коло
Результати
Висновок
В цій лабораторній роботі я оцінила якість генераторів ПВЧ. Оцінка якості методів генерації ПВЧ за методом гістограм показала, що найгіршим є метод середніх квадратів, а найефективнішим – мультиплікативний метод, а стандартний генератор виявився кращим за вихор Мерсенна.
При перевірці за непрямими ознаками також найкращим виявився Вихор Мерсенна, а найгіршим – метод середніх квадратів.
По довжині періоду найкращим виявився вихор Мерсенна, а найбільший період серед інших методів мав мультиплікативний метод.