libdspl-2.0
Digital Signal Processing Algorithm Library
randgen.c
1 /*
2 * Copyright (c) 2015-2019 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of libdspl-2.0.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * DSPL is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #include <stdlib.h>
23 #include <math.h>
24 #include <time.h>
25 
26 #include "dspl.h"
27 #include "dspl_internal.h"
28 #include "mt19937.h"
29 
30 
31 
32 #ifdef DOXYGEN_ENGLISH
33 
75 #endif
76 #ifdef DOXYGEN_RUSSIAN
77 
120 #endif
121 int DSPL_API random_init(random_t* prnd, int type, void* seed)
122 {
123  srand(time(NULL));
124 
125  if(!prnd)
126  return RES_OK;
127 
128  switch(type)
129  {
130  case RAND_TYPE_MRG32K3A:
131  /* MRG32k3a init */
132  prnd->mrg32k3a_x[0] = prnd->mrg32k3a_x[1] = 1.0;
133  prnd->mrg32k3a_y[0] = prnd->mrg32k3a_y[1] =
134  prnd->mrg32k3a_y[2] = 1.0;
135  if(seed)
136  prnd->mrg32k3a_x[2] = *((double*)seed);
137  else
138  prnd->mrg32k3a_x[2] = (double) rand() * rand();
139  break;
140  case RAND_TYPE_MT19937:
141  if(seed)
142  mt19937_init_genrand64(*((unsigned long long*)seed), prnd);
143  else
144  mt19937_init_genrand64((unsigned long long)rand()*rand(), prnd);
145  break;
146  default:
147  return ERROR_RAND_TYPE;
148  }
149  prnd->type = type;
150 
151  return RES_OK;
152 }
153 
154 
155 
156 
157 #ifdef DOXYGEN_ENGLISH
158 
199 #endif
200 #ifdef DOXYGEN_RUSSIAN
201 
243 #endif
244 int DSPL_API randb(double* x, int n, random_t* prnd)
245 {
246  double z[RAND_BUFSIZE];
247  int i, cnt, err;
248  if(!x)
249  return ERROR_PTR;
250  if(n < 1)
251  return ERROR_SIZE;
252  cnt = 0;
253  while(cnt < n)
254  {
255  i = cnt % RAND_BUFSIZE;
256  if(!i)
257  {
258  err = randu(z, RAND_BUFSIZE, prnd);
259  if(err != RES_OK)
260  return err;
261  }
262  x[cnt] = z[i] > 0.5 ? 1.0 : 0.0;
263  cnt++;
264  }
265  return RES_OK;
266 }
267 
268 
269 
270 
271 #ifdef DOXYGEN_ENGLISH
272 
313 #endif
314 #ifdef DOXYGEN_RUSSIAN
315 
357 #endif
358 int DSPL_API randb2(double* x, int n, random_t* prnd)
359 {
360  double z[RAND_BUFSIZE];
361  int i, cnt, err;
362  if(!x)
363  return ERROR_PTR;
364  if(n < 1)
365  return ERROR_SIZE;
366  cnt = 0;
367  while(cnt < n)
368  {
369  i = cnt % RAND_BUFSIZE;
370  if(!i)
371  {
372  err = randu(z, RAND_BUFSIZE, prnd);
373  if(err != RES_OK)
374  return err;
375  }
376  x[cnt] = z[i] > 0.5 ? 1.0 : -1.0;
377  cnt++;
378  }
379  return RES_OK;
380 }
381 
382 
383 
384 
385 #ifdef DOXYGEN_ENGLISH
386 
387 #endif
388 #ifdef DOXYGEN_RUSSIAN
389 
390 #endif
391 int randu_mrg32k3a (double* u, int n, random_t* prnd)
392 {
393 
394  if(!u || !prnd)
395  return ERROR_PTR;
396  if(n < 1)
397  return ERROR_SIZE;
398 
399  long z;
400  double xn, yn, *x, *y;
401  int k;
402 
403  x = prnd->mrg32k3a_x;
404  y = prnd->mrg32k3a_y;
405  for(k = 0; k < n; k++)
406  {
407  /* Component x[n] */
408  xn = MRG32K3A_A12 * x[1] - MRG32K3A_A13 * x[2];
409 
410  z = (long)(xn / MRG32K3A_M1);
411  xn -= (double)z * MRG32K3A_M1;
412  if (xn < 0.0)
413  xn += MRG32K3A_M1;
414 
415  x[2] = x[1];
416  x[1] = x[0];
417  x[0] = xn;
418 
419  /* Component y[n] */
420  yn = MRG32K3A_A21 * y[0] - MRG32K3A_A23 * y[2];
421  z = (long)(yn / MRG32K3A_M2);
422  yn -= (double)z * MRG32K3A_M2;
423  if (yn < 0.0)
424  yn += MRG32K3A_M2;
425 
426  y[2] = y[1];
427  y[1] = y[0];
428  y[0] = yn;
429 
430  /* Combination */
431  u[k] = (xn <= yn) ? ((xn - yn + MRG32K3A_M1) * MRG32K3A_NORM):
432  (xn - yn) * MRG32K3A_NORM;
433  }
434  return RES_OK;
435 }
436 
437 
438 
439 
440 
441 
442 
443 
444 #ifdef DOXYGEN_ENGLISH
445 
446 #endif
447 #ifdef DOXYGEN_RUSSIAN
448 
497 #endif
498 int DSPL_API randi(int* x, int n, int start, int stop, random_t* prnd)
499 {
500  double z[RAND_BUFSIZE];
501  double dx;
502  int i, cnt, err;
503  if(!x)
504  return ERROR_PTR;
505  if(n < 1)
506  return ERROR_SIZE;
507 
508  dx = (double)stop - (double)start;
509  cnt = 0;
510  while(cnt < n)
511  {
512  i = cnt % RAND_BUFSIZE;
513  if(!i)
514  {
515  err = randu(z, RAND_BUFSIZE, prnd);
516  if(err != RES_OK)
517  return err;
518  }
519  x[cnt] = start + (int)round(z[i] * dx);
520  cnt++;
521  }
522  return RES_OK;
523 }
524 
525 
526 
527 
528 
529 
530 
531 #ifdef DOXYGEN_ENGLISH
532 
533 #endif
534 #ifdef DOXYGEN_RUSSIAN
535 
588 #endif
589 int DSPL_API randn(double* x, int n, double mu, double sigma, random_t* prnd)
590 {
591  int k, m;
592  double x1[RAND_BUFSIZE], x2[RAND_BUFSIZE];
593  int res;
594  if(!x)
595  return ERROR_PTR;
596 
597  if(n<1)
598  return ERROR_SIZE;
599 
600  if(sigma < 0.0)
601  return ERROR_RAND_SIGMA;
602 
603  k=0;
604  while(k < n)
605  {
606  if((res = randu(x1, RAND_BUFSIZE, prnd)) != RES_OK)
607  goto exit_label;
608  if((res = randu(x2, RAND_BUFSIZE, prnd)) != RES_OK)
609  goto exit_label;
610  m = 0;
611  while(k < n && m < RAND_BUFSIZE)
612  {
613  if(x1[m] != 0.0)
614  {
615  x[k] = sqrt(-2.0*log(x1[m]))*cos(M_2PI*x2[m])*sigma + mu;
616  k++;
617  m++;
618  }
619  }
620  }
621 
622  res = RES_OK;
623 exit_label:
624  return res;
625 }
626 
627 
628 
629 
630 #ifdef DOXYGEN_ENGLISH
631 
632 #endif
633 #ifdef DOXYGEN_RUSSIAN
634 
635 #endif
636 int DSPL_API randn_cmplx(complex_t* x, int n, complex_t* mu,
637  double sigma, random_t* prnd)
638 {
639  int err, i;
640 
641  err = randn((double*)x, 2*n, 0.0, sigma / M_SQRT2, prnd);
642  if(err!= RES_OK)
643  return err;
644  if(mu)
645  {
646  for(i = 0; i < n; i++)
647  {
648  RE(x[i]) += RE(mu[0]);
649  IM(x[i]) += IM(mu[0]);
650  }
651  }
652  return err;
653 }
654 
655 
656 
657 
658 
659 
660 #ifdef DOXYGEN_ENGLISH
661 
662 #endif
663 #ifdef DOXYGEN_RUSSIAN
664 
726 #endif
727 int DSPL_API randu(double* x, int n, random_t* prnd)
728 {
729  int i;
730 
731  if(!x)
732  return ERROR_PTR;
733  if(n < 0)
734  return ERROR_SIZE;
735 
736  if(prnd)
737  {
738  switch(prnd->type)
739  {
740  case RAND_TYPE_MRG32K3A:
741  return randu_mrg32k3a(x, n, prnd);
742  case RAND_TYPE_MT19937:
743  for(i = 0; i < n; i++)
744  x[i] = mt19937_genrand64_real1(prnd);
745  return RES_OK;
746  default:
747  return ERROR_RAND_TYPE;
748  }
749  }
750  else
751  {
752  if(!x)
753  return ERROR_PTR;
754  if(n<1)
755  return ERROR_SIZE;
756  for(i = 0; i < n; i++)
757  x[i] = (double)rand()/RAND_MAX;
758  }
759 
760  return RES_OK;
761 }
762 
random_t
Definition: dspl.h:350
RE
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:420
ERROR_PTR
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:610
ERROR_SIZE
#define ERROR_SIZE
Error array size. This error occurs when in addition to the pointer the wrong input is passed to the ...
Definition: dspl.h:618
random_t::mrg32k3a_y
double mrg32k3a_y[3]
Definition: dspl.h:354
random_t::type
int type
Definition: dspl.h:360
ERROR_RAND_TYPE
#define ERROR_RAND_TYPE
Unknown pseudorandom number genration algorithm. The following algorithms are used in the library:
Definition: dspl.h:614
random_init
int DSPL_API random_init(random_t *prnd, int type, void *seed)
Pseudorandom numbers generators initialization.
Definition: randgen.c:121
randb
int DSPL_API randb(double *x, int n, random_t *prnd)
Binary unipolar [0, 1] pseudorandom vector.
Definition: randgen.c:244
complex_t
double complex_t[2]
Complex data type.
Definition: dspl.h:86
RES_OK
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:558
random_t::mrg32k3a_x
double mrg32k3a_x[3]
Definition: dspl.h:353
IM
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:478
ERROR_RAND_SIGMA
#define ERROR_RAND_SIGMA
The standard deviation is incorrect normal distribution of a random variable. The standard deviation ...
Definition: dspl.h:613
randb2
int DSPL_API randb2(double *x, int n, random_t *prnd)
Binary bipolar [-1, 1] pseudorandom vector.
Definition: randgen.c:358
Generated on Mon Dec 14 2020 14:02:36 for libdspl-2.0 by   doxygen 1.8.20
libdspl-2.0 documentation is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.