summaryrefslogtreecommitdiff
path: root/packages/base
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-12-01 17:11:55 +0100
committerAlberto Ruiz <aruiz@um.es>2014-12-01 17:11:55 +0100
commit27c9e9b804d7e55f5ae180040f76d58116a85b08 (patch)
tree4d21cedb1e8e72975ce73ad96fa6f50006226fc1 /packages/base
parentbb7dce54060fe60bdb04af9f910680cc460d4d5b (diff)
remove static state in gaussrand
Diffstat (limited to 'packages/base')
-rw-r--r--packages/base/src/C/vector-aux.c16
1 files changed, 10 insertions, 6 deletions
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index d5b6d4d..f8feb37 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -708,13 +708,13 @@ inline double urandom(struct drand48_data * buffer) {
708 708
709 709
710// http://c-faq.com/lib/gaussian.html 710// http://c-faq.com/lib/gaussian.html
711double gaussrand(struct drand48_data * buffer) 711double gaussrand(struct drand48_data *buffer,
712 int *phase, double *pV1, double *pV2, double *pS)
712{ 713{
713 static double V1, V2, S; 714 double V1=*pV1, V2=*pV2, S=*pS;
714 static int phase = 0;
715 double X; 715 double X;
716 716
717 if(phase == 0) { 717 if(*phase == 0) {
718 do { 718 do {
719 double U1 = urandom(buffer); 719 double U1 = urandom(buffer);
720 double U2 = urandom(buffer); 720 double U2 = urandom(buffer);
@@ -728,7 +728,8 @@ double gaussrand(struct drand48_data * buffer)
728 } else 728 } else
729 X = V2 * sqrt(-2 * log(S) / S); 729 X = V2 * sqrt(-2 * log(S) / S);
730 730
731 phase = 1 - phase; 731 *phase = 1 - *phase;
732 *pV1=V1; *pV2=V2; *pS=S;
732 733
733 return X; 734 return X;
734} 735}
@@ -736,6 +737,9 @@ double gaussrand(struct drand48_data * buffer)
736int random_vector(unsigned int seed, int code, DVEC(r)) { 737int random_vector(unsigned int seed, int code, DVEC(r)) {
737 struct drand48_data buffer; 738 struct drand48_data buffer;
738 srand48_r(seed,&buffer); 739 srand48_r(seed,&buffer);
740 int phase = 0;
741 double V1,V2,S;
742
739 int k; 743 int k;
740 switch (code) { 744 switch (code) {
741 case 0: { // uniform 745 case 0: { // uniform
@@ -746,7 +750,7 @@ int random_vector(unsigned int seed, int code, DVEC(r)) {
746 } 750 }
747 case 1: { // gaussian 751 case 1: { // gaussian
748 for (k=0; k<rn; k++) { 752 for (k=0; k<rn; k++) {
749 rp[k] = gaussrand(&buffer); 753 rp[k] = gaussrand(&buffer,&phase,&V1,&V2,&S);
750 } 754 }
751 OK 755 OK
752 } 756 }