merge master

3 files changed, 80 insertions, 5 deletions

diff --git a/packages/base/THANKS.md b/packages/base/THANKS.md
index 90e6199..c220f16 100644
--- a/packages/base/THANKS.md
+++ b/packages/base/THANKS.md
@@ -177,3 +177,5 @@ module reorganization, monadic mapVectorM, and many other improvements.
 
 - "yongqli" reported the bug in randomVector (rand() is not thread-safe).
 
+- "yongqli" and Kiwamu Ishikura reported that drand48_r() is not portable.
+
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index 2900c42..b3559a9 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.16.1.0
+Version:            0.16.1.2
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -118,6 +118,7 @@ library
         if arch(i386)
             cc-options: -arch i386
         frameworks: Accelerate
+        cpp-options: -DOSX
 
     if os(freebsd)
        extra-lib-dirs: /usr/local/lib
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index d01e788..0f1b3c7 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -700,7 +700,10 @@ int saveMatrix(char * file, char * format, KDMAT(a)){
 
 ////////////////////////////////////////////////////////////////////////////////
 
-// http://c-faq.com/lib/gaussian.html
+#ifdef OSX
+
+#pragma message "randomVector is not thread-safe in OSX"
+
 double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
 {
         double V1=*pV1, V2=*pV2, S=*pS;
@@ -730,10 +733,8 @@ 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;
+    
     int k;
-
-    srandom(seed);
-
     switch (code) {
       case 0: { // uniform
         for (k=0; k<rn; k++) {
@@ -752,6 +753,77 @@ int random_vector(unsigned int seed, int code, DVEC(r)) {
     }
 }
 
+#else
+
+inline double urandom(struct random_data * buffer) {
+    int32_t res;
+    random_r(buffer,&res);
+    return (double)res/RAND_MAX;
+}
+
+
+// http://c-faq.com/lib/gaussian.html
+double gaussrand(struct random_data *buffer,
+                 int *phase, double *pV1, double *pV2, double *pS)
+{
+        double V1=*pV1, V2=*pV2, S=*pS;
+        double X;
+
+        if(*phase == 0) {
+                do {
+            double U1 = (double)random() / (double)RAND_MAX;
+                        double U2 = (double)random() / (double)RAND_MAX;
+
+                        V1 = 2 * U1 - 1;
+                        V2 = 2 * U2 - 1;
+                        S = V1 * V1 + V2 * V2;
+                        } while(S >= 1 || S == 0);
+
+                X = V1 * sqrt(-2 * log(S) / S);
+        } else
+                X = V2 * sqrt(-2 * log(S) / S);
+
+        *phase = 1 - *phase;
+    *pV1=V1; *pV2=V2; *pS=S;
+
+        return X;
+
+}
+
+int random_vector(unsigned int seed, int code, DVEC(r)) {
+    struct random_data buffer;
+    char   random_state[128];
+    memset(&buffer, 0, sizeof(struct random_data));
+    memset(random_state, 0, sizeof(random_state));
+    
+    initstate_r(seed,random_state,sizeof(random_state),&buffer);
+    // setstate_r(random_state,&buffer);
+    // srandom_r(seed,&buffer);
+
+    int phase = 0;
+    double V1,V2,S;
+
+    int k;
+    switch (code) {
+      case 0: { // uniform
+        for (k=0; k<rn; k++) {
+            rp[k] = urandom(&buffer);
+        }
+        OK
+      }
+      case 1: { // gaussian
+        for (k=0; k<rn; k++) {
+            rp[k] = gaussrand(&buffer,&phase,&V1,&V2,&S);
+        }
+        OK
+      }
+
+      default: ERROR(BAD_CODE);
+    }
+}
+
+#endif
+
 ////////////////////////////////////////////////////////////////////////////////
 
 int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {