From 163dd4d4e823d74d6175664abec75d8819930b3e Mon Sep 17 00:00:00 2001 From: Dsplib Date: Sun, 31 Mar 2019 20:08:11 +0300 Subject: [PATCH] added mrg32k3a random algorithm --- dspl/src/dspl_internal.h | 12 ++++++- dspl/src/randgen.c | 74 +++++++++++++++++++++++++++++++++++---- include/dspl.c | 2 ++ include/dspl.h | 21 +++++++++-- release/include/dspl.c | 2 ++ release/include/dspl.h | 21 +++++++++-- verif/src/readbin_verif.c | 7 ++-- 7 files changed, 125 insertions(+), 14 deletions(-) diff --git a/dspl/src/dspl_internal.h b/dspl/src/dspl_internal.h index e87085a..a24e62c 100644 --- a/dspl/src/dspl_internal.h +++ b/dspl/src/dspl_internal.h @@ -89,8 +89,18 @@ int win_rect (double *w, int n); int fir_linphase_lpf(int ord, double wp, int wintype, double winparam, double* h); -#define MATRIX_SINGULAR_THRESHOLD 1E-14 +#define MATRIX_SINGULAR_THRESHOLD 1E-14 +/* random MRG32K3A algorithm constants */ +#define MRG32K3A_NORM 2.328306549295728e-10 +#define MRG32K3A_M1 4294967087.0 +#define MRG32K3A_M2 4294944443.0 +#define MRG32K3A_A12 1403580.0 +#define MRG32K3A_A13 810728.0 +#define MRG32K3A_A21 527612.0 +#define MRG32K3A_A23 1370589.0 + + #endif diff --git a/dspl/src/randgen.c b/dspl/src/randgen.c index be00189..4be4c37 100644 --- a/dspl/src/randgen.c +++ b/dspl/src/randgen.c @@ -29,12 +29,72 @@ + +void DSPL_API random_init(random_t* prnd) +{ + srand(time(NULL)); + // MRG32k3a init + prnd->mrg32k3a_x[0] = prnd->mrg32k3a_x[1] = 1.0; + prnd->mrg32k3a_y[0] = prnd->mrg32k3a_y[1] = prnd->mrg32k3a_y[2] = 1.0; + prnd->mrg32k3a_x[2] = rand(); +} + +int randu_mrg32k3a (double* u, int n, random_t* prnd) +{ + + if(!u || !prnd) + return ERROR_PTR; + if(n < 1) + return ERROR_SIZE; + + long z; + double xn, yn, *x, *y; + int k; + + x = prnd->mrg32k3a_x; + y = prnd->mrg32k3a_y; + for(k = 0; k < n; k++) + { + /* Component x[n] */ + xn = MRG32K3A_A12 * x[1] - MRG32K3A_A13 * x[2]; + + z = (long)(xn / MRG32K3A_M1); + xn -= (double)z * MRG32K3A_M1; + if (xn < 0.0) + xn += MRG32K3A_M1; + + x[2] = x[1]; + x[1] = x[0]; + x[0] = xn; + + /* Component y[n] */ + yn = MRG32K3A_A21 * y[0] - MRG32K3A_A23 * y[2]; + z = (long)(yn / MRG32K3A_M2); + yn -= (double)z * MRG32K3A_M2; + if (yn < 0.0) + yn += MRG32K3A_M2; + + y[2] = y[1]; + y[1] = y[0]; + y[0] = yn; + + /* Combination */ + u[k] = (xn <= yn) ? ((xn - yn + MRG32K3A_M1) * MRG32K3A_NORM): + (xn - yn) * MRG32K3A_NORM; + } + return RES_OK; +} + + + + /****************************************************************************** Uniform random numbers generator *******************************************************************************/ -int DSPL_API randu(double* x, int n) +int DSPL_API randu(double* x, int n, random_t* prnd) { - int k,m; + + /*int k,m; unsigned int x1[4], x2[4], y; if(!x) @@ -61,8 +121,8 @@ int DSPL_API randu(double* x, int n) x[k] = (double)y/DSPL_RAND_MOD_X1; } - - return RES_OK; + */ + return randu_mrg32k3a(x, n, prnd); } @@ -72,7 +132,7 @@ int DSPL_API randu(double* x, int n) /******************************************************************************* Gaussian random numbers generator *******************************************************************************/ -int DSPL_API randn(double* x, int n, double mu, double sigma) +int DSPL_API randn(double* x, int n, double mu, double sigma, random_t* prnd) { int k, m; double x1[128], x2[128]; @@ -89,11 +149,11 @@ int DSPL_API randn(double* x, int n, double mu, double sigma) k=0; while(k < n) { - res = randu(x1, 128); + res = randu(x1, 128, prnd); if(res != RES_OK) goto exit_label; - res = randu(x2, 128); + res = randu(x2, 128, prnd); if(res != RES_OK) goto exit_label; m = 0 ; diff --git a/include/dspl.c b/include/dspl.c index 1c37b91..91eb97c 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -116,6 +116,7 @@ p_poly_z2a_cmplx poly_z2a_cmplx ; p_polyval polyval ; p_polyval_cmplx polyval_cmplx ; p_randn randn ; +p_random_init random_init ; p_randu randu ; p_ratcompos ratcompos ; p_re2cmplx re2cmplx ; @@ -271,6 +272,7 @@ void* dspl_load() LOAD_FUNC(polyval); LOAD_FUNC(polyval_cmplx); LOAD_FUNC(randn); + LOAD_FUNC(random_init); LOAD_FUNC(randu); LOAD_FUNC(ratcompos); LOAD_FUNC(re2cmplx); diff --git a/include/dspl.h b/include/dspl.h index 670dfd1..e7bad81 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -42,6 +42,7 @@ typedef double complex_t[2]; + typedef struct { complex_t* w; @@ -61,6 +62,18 @@ typedef struct } matrix_t; + +typedef struct +{ + double mrg32k3a_seed; + double mrg32k3a_x[3]; + double mrg32k3a_y[3]; + +}random_t; + + + + #define RE(x) (x[0]) #define IM(x) (x[1]) @@ -672,10 +685,14 @@ DECLARE_FUNC(int, polyval_cmplx, complex_t* DECLARE_FUNC(int, randn, double* COMMA int COMMA double - COMMA double); + COMMA double + COMMA random_t* prnd); +//------------------------------------------------------------------------------ +DECLARE_FUNC(void, random_init, random_t* prnd); //------------------------------------------------------------------------------ DECLARE_FUNC(int, randu, double* - COMMA int); + COMMA int + COMMA random_t* prnd); //------------------------------------------------------------------------------ DECLARE_FUNC(int, ratcompos, double* b COMMA double* a diff --git a/release/include/dspl.c b/release/include/dspl.c index 1c37b91..91eb97c 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -116,6 +116,7 @@ p_poly_z2a_cmplx poly_z2a_cmplx ; p_polyval polyval ; p_polyval_cmplx polyval_cmplx ; p_randn randn ; +p_random_init random_init ; p_randu randu ; p_ratcompos ratcompos ; p_re2cmplx re2cmplx ; @@ -271,6 +272,7 @@ void* dspl_load() LOAD_FUNC(polyval); LOAD_FUNC(polyval_cmplx); LOAD_FUNC(randn); + LOAD_FUNC(random_init); LOAD_FUNC(randu); LOAD_FUNC(ratcompos); LOAD_FUNC(re2cmplx); diff --git a/release/include/dspl.h b/release/include/dspl.h index 670dfd1..e7bad81 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -42,6 +42,7 @@ typedef double complex_t[2]; + typedef struct { complex_t* w; @@ -61,6 +62,18 @@ typedef struct } matrix_t; + +typedef struct +{ + double mrg32k3a_seed; + double mrg32k3a_x[3]; + double mrg32k3a_y[3]; + +}random_t; + + + + #define RE(x) (x[0]) #define IM(x) (x[1]) @@ -672,10 +685,14 @@ DECLARE_FUNC(int, polyval_cmplx, complex_t* DECLARE_FUNC(int, randn, double* COMMA int COMMA double - COMMA double); + COMMA double + COMMA random_t* prnd); +//------------------------------------------------------------------------------ +DECLARE_FUNC(void, random_init, random_t* prnd); //------------------------------------------------------------------------------ DECLARE_FUNC(int, randu, double* - COMMA int); + COMMA int + COMMA random_t* prnd); //------------------------------------------------------------------------------ DECLARE_FUNC(int, ratcompos, double* b COMMA double* a diff --git a/verif/src/readbin_verif.c b/verif/src/readbin_verif.c index 35670e8..621dbb3 100644 --- a/verif/src/readbin_verif.c +++ b/verif/src/readbin_verif.c @@ -17,8 +17,11 @@ int main() xr = (double*) malloc(N * sizeof(double)); xc = (complex_t*) malloc(N * sizeof(complex_t)); - randn(xr, N, 0, 1.0); - randn((double*)xc, 2*N, 0, 1.0); + random_t rnd; + random_init(&rnd); + + randn(xr, N, 0, 1.0,&rnd); + randn((double*)xc, 2*N, 0, 1.0,&rnd); writebin(xr, N, DAT_DOUBLE, "dat/in.dat");