kopia lustrzana https://github.com/Dsplib/libdspl-2.0
added mrg32k3a random algorithm
Dsplib
rodzic
516fa73ace
commit
163dd4d4e8
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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 ;
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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");
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue