I tried to simulate normal random number. From online source recommend, xoshiro256** faster than rand() in c++. My idea is to first simulate uniform RN and transform into normal rn
#define _USE_MATH_DEFINES
#include <iostream>
#include <stdint.h>
#include <cmath>
#include <chrono>
double GetOneGaussianByBoxMuller()
{
double result;
double x;
double y;
double sizeSquared;
do
{
x = 2.0*rand()/static_cast<double>(RAND_MAX)-1;
y = 2.0*rand()/static_cast<double>(RAND_MAX)-1;
sizeSquared = x*x + y*y;
}
while( sizeSquared >= 1.0 || sizeSquared ==0);
result = x*sqrt(-2*log(sizeSquared)/sizeSquared);
return result;
};
uint64_t rol64(uint64_t x, int k)
{
return (x << k) | (x >> (64 - k));
};
struct xoshiro256ss_state {
uint64_t s[4] = {5, 5, 2, 1};
};
uint64_t xoshiro256ss(struct xoshiro256ss_state *state)
{
uint64_t *s = state->s;
uint64_t const result = rol64(s[1] * 5, 7) * 9;
uint64_t const t = s[1] << 17;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = rol64(s[3], 45);
for (int i=0; i<4; i++){
std::cout<<s[i]<<std::endl;
}
return result;
};
int main(){
const int m = 1000;
double rand_max = pow(2,64);
xoshiro256ss_state seed;
double result;
double x;
double y;
double sizeSquared;
auto t1 = std::chrono::high_resolution_clock::now();
for (int i = 0; i < m ; i++){
do
{
x = 2.0* xoshiro256ss(&seed) / rand_max-1;
y = 2.0* xoshiro256ss(&seed) / rand_max-1;
sizeSquared = x*x + y*y;
//std::cout<< x <<std::endl;
//std::cout<< y <<std::endl;
}
while( sizeSquared >= 1.0 || sizeSquared ==0);
std::cout<< x*sqrt(-2*log(sizeSquared)/sizeSquared) <<std::endl;
}
auto t2 = std::chrono::high_resolution_clock::now();
auto t3 = std::chrono::high_resolution_clock::now();
for (int i = 0; i < m ; i++){
std::cout<< GetOneGaussianByBoxMuller() <<std::endl;
}
auto t4 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> ms_double = t2 - t1;
std::chrono::duration<double, std::milli> ms_double1 = t4 - t3;
std::cout << ms_double.count() << "ms";
std::cout << ms_double1.count() << "ms";
}
The bitwise operation in function uint64_t xoshiro256ss(struct xoshiro256ss_state *state)
seems not working when it is too big. Also, it seems like my code about xoshiro256** is wrong as well. Can anyone help me with this.
Thanks
Aucun commentaire:
Enregistrer un commentaire