diff --git a/THANKS.md b/THANKS.md
--- a/THANKS.md
+++ b/THANKS.md
@@ -175,7 +175,7 @@
 
 - John Shahbazian added support for openBLAS.
 
-- "yongqli" reported the bug in randomVector (rand() is not thread-safe).
+- "yongqli" reported the bug in randomVector (rand() is not thread-safe and drand48_r() is not portable).
 
-- "yongqli" and Kiwamu Ishikura reported that drand48_r() is not portable.
+- Kiwamu Ishikura improved randomVector for OSX
 
diff --git a/hmatrix.cabal b/hmatrix.cabal
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.16.1.2
+Version:            0.16.1.3
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -118,7 +118,6 @@
         if arch(i386)
             cc-options: -arch i386
         frameworks: Accelerate
-        cpp-options: -DOSX
 
     if os(freebsd)
        extra-lib-dirs: /usr/local/lib
diff --git a/src/C/vector-aux.c b/src/C/vector-aux.c
--- a/src/C/vector-aux.c
+++ b/src/C/vector-aux.c
@@ -700,10 +700,15 @@
 
 ////////////////////////////////////////////////////////////////////////////////
 
-#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 @@
 
 	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 @@
 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 @@
 
 	if(*phase == 0) {
 		do {
-			double U1 = urandom(buffer);
+            double U1 = urandom(buffer);
 			double U2 = urandom(buffer);
 
 			V1 = 2 * U1 - 1;
@@ -787,6 +794,7 @@
     *pV1=V1; *pV2=V2; *pS=S;
 
 	return X;
+
 }
 
 int random_vector(unsigned int seed, int code, DVEC(r)) {
@@ -801,7 +809,7 @@
 
     int phase = 0;
     double V1,V2,S;
-    
+
     int k;
     switch (code) {
       case 0: { // uniform