dimanche 18 avril 2021

random number from xoshiro256** to box-muller

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