Merge pull request #101 from i-kiwamu/use-random

use random() instead of drand48_r

1 files changed, 14 insertions, 6 deletions

diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index 946d069..dda47cb 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -700,10 +700,15 @@ int saveMatrix(char * file, char * format, KDMAT(a)){
 
 ////////////////////////////////////////////////////////////////////////////////
 
-#ifdef OSX
+#ifdef __APPLE__
 
 #pragma message "randomVector is not thread-safe in OSX"
 
+inline double urandom() {
+    const long max_random = 2147483647; // 2**31 - 1
+    return (double)random() / (double)max_random;
+}
+
 double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
 {
         double V1=*pV1, V2=*pV2, S=*pS;
@@ -711,8 +716,8 @@ double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
 
         if(*phase == 0) {
                 do {
-            double U1 = (double)random() / (double)RAND_MAX;
-                        double U2 = (double)random() / (double)RAND_MAX;
+            double U1 = urandom();
+                        double U2 = urandom();
 
                         V1 = 2 * U1 - 1;
                         V2 = 2 * U2 - 1;
@@ -733,12 +738,14 @@ double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
 int random_vector(unsigned int seed, int code, DVEC(r)) {
     int phase = 0;
     double V1,V2,S;
+
+    srandom(seed);
     
     int k;
     switch (code) {
       case 0: { // uniform
         for (k=0; k<rn; k++) {
-            rp[k] = (double)random() / (double)RAND_MAX;
+            rp[k] = urandom();
         }
         OK
       }
@@ -771,7 +778,7 @@ double gaussrand(struct random_data *buffer,
 
         if(*phase == 0) {
                 do {
-                        double U1 = urandom(buffer);
+            double U1 = urandom(buffer);
                         double U2 = urandom(buffer);
 
                         V1 = 2 * U1 - 1;
@@ -787,6 +794,7 @@ double gaussrand(struct random_data *buffer,
     *pV1=V1; *pV2=V2; *pS=S;
 
         return X;
+
 }
 
 int random_vector(unsigned int seed, int code, DVEC(r)) {
@@ -801,7 +809,7 @@ int random_vector(unsigned int seed, int code, DVEC(r)) {
 
     int phase = 0;
     double V1,V2,S;
-    
+
     int k;
     switch (code) {
       case 0: { // uniform