summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--packages/base/src/C/vector-aux.c17
1 files changed, 12 insertions, 5 deletions
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index 0f1b3c7..1c3fe59 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -704,6 +704,11 @@ int saveMatrix(char * file, char * format, KDMAT(a)){
704 704
705#pragma message "randomVector is not thread-safe in OSX" 705#pragma message "randomVector is not thread-safe in OSX"
706 706
707inline double urandom() {
708 const long max_random = 2147483647 // 2**31 - 1
709 return (double)random() / (double)max_random;
710}
711
707double gaussrand(int *phase, double *pV1, double *pV2, double *pS) 712double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
708{ 713{
709 double V1=*pV1, V2=*pV2, S=*pS; 714 double V1=*pV1, V2=*pV2, S=*pS;
@@ -711,8 +716,8 @@ double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
711 716
712 if(*phase == 0) { 717 if(*phase == 0) {
713 do { 718 do {
714 double U1 = (double)random() / (double)RAND_MAX; 719 double U1 = urandom();
715 double U2 = (double)random() / (double)RAND_MAX; 720 double U2 = urandom();
716 721
717 V1 = 2 * U1 - 1; 722 V1 = 2 * U1 - 1;
718 V2 = 2 * U2 - 1; 723 V2 = 2 * U2 - 1;
@@ -733,12 +738,14 @@ double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
733int random_vector(unsigned int seed, int code, DVEC(r)) { 738int random_vector(unsigned int seed, int code, DVEC(r)) {
734 int phase = 0; 739 int phase = 0;
735 double V1,V2,S; 740 double V1,V2,S;
741
742 srandom(seed);
736 743
737 int k; 744 int k;
738 switch (code) { 745 switch (code) {
739 case 0: { // uniform 746 case 0: { // uniform
740 for (k=0; k<rn; k++) { 747 for (k=0; k<rn; k++) {
741 rp[k] = (double)random() / (double)RAND_MAX; 748 rp[k] = urandom();
742 } 749 }
743 OK 750 OK
744 } 751 }
@@ -771,8 +778,8 @@ double gaussrand(struct random_data *buffer,
771 778
772 if(*phase == 0) { 779 if(*phase == 0) {
773 do { 780 do {
774 double U1 = (double)random() / (double)RAND_MAX; 781 double U1 = urandom(buffer);
775 double U2 = (double)random() / (double)RAND_MAX; 782 double U2 = urandom(buffer);
776 783
777 V1 = 2 * U1 - 1; 784 V1 = 2 * U1 - 1;
778 V2 = 2 * U2 - 1; 785 V2 = 2 * U2 - 1;