diff --git a/CUtils/c_sources/mt19937ar.c b/CUtils/c_sources/mt19937ar.c
index 3c99f71e619566f7a02e8f07cdc8276b44d19ee5..af5e5f117d55424a6cd27a683b2b921a88f6bf0f 100644
--- a/CUtils/c_sources/mt19937ar.c
+++ b/CUtils/c_sources/mt19937ar.c
@@ -42,6 +42,7 @@
 */
 
 #include <stdio.h>
+#include "debug.h"
 #include "mt19937ar.h"
 
 /* Period parameters */  
@@ -51,7 +52,7 @@
 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
 
-static mt19937ar_t local = { {0, }, /* the array for the state vector  */
+static __thread mt19937ar_t local = { {0, }, /* the array for the state vector  */
 			     N+1, /* mti==N+1 means mt[N] is not initialized */
 };
 
@@ -60,6 +61,7 @@ void init_genrand_mt(mt19937ar_t* p, unsigned long s)
 {
     int mti;
     unsigned long *mt=&p->mt[0];
+    debug("init init_genrand_mt(%li)", s);
     mt[0]= s & 0xffffffffUL;
     for (mti=1; mti<N; mti++) {
         mt[mti] = 
@@ -87,6 +89,7 @@ void init_by_array_mt(mt19937ar_t* p, unsigned long init_key[], int key_length)
 {
     int i, j, k;
     unsigned long *mt=&p->mt[0];
+    debug("init init_by_array_mt");
     init_genrand_mt(p, 19650218UL);
     i=1; j=0;
     k = (N>key_length ? N : key_length);
@@ -163,6 +166,11 @@ long genrand_int31(void)
     return (long)(genrand_int32()>>1);
 }
 
+long genrand_int31_mt(mt19937ar_t* p)
+{
+    return (long)(genrand_int32_mt(p)>>1);
+}
+
 /* generates a random number on [0,1]-real-interval */
 double genrand_real1(void)
 {
@@ -170,6 +178,12 @@ double genrand_real1(void)
     /* divided by 2^32-1 */ 
 }
 
+double genrand_real1_mt(mt19937ar_t* p)
+{
+    return genrand_int32_mt(p)*(1.0/4294967295.0); 
+    /* divided by 2^32-1 */ 
+}
+
 /* generates a random number on [0,1)-real-interval */
 double genrand_real2(void)
 {
@@ -177,6 +191,12 @@ double genrand_real2(void)
     /* divided by 2^32 */
 }
 
+double genrand_real2_mt(mt19937ar_t* p)
+{
+    return genrand_int32_mt(p)*(1.0/4294967296.0); 
+    /* divided by 2^32 */
+}
+
 /* generates a random number on (0,1)-real-interval */
 double genrand_real3(void)
 {
@@ -184,12 +204,24 @@ double genrand_real3(void)
     /* divided by 2^32 */
 }
 
+double genrand_real3_mt(mt19937ar_t* p)
+{
+    return (((double)genrand_int32_mt(p)) + 0.5)*(1.0/4294967296.0); 
+    /* divided by 2^32 */
+}
+
 /* generates a random number on [0,1) with 53-bit resolution*/
 double genrand_res53(void) 
 { 
     unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; 
     return(a*67108864.0+b)*(1.0/9007199254740992.0); 
 } 
+
+double genrand_res53_mt(mt19937ar_t* p) 
+{ 
+    unsigned long a=genrand_int32_mt(p)>>5, b=genrand_int32_mt(p)>>6; 
+    return(a*67108864.0+b)*(1.0/9007199254740992.0); 
+} 
 /* These real versions are due to Isaku Wada, 2002/01/09 added */
 
 #if 0
diff --git a/CUtils/c_sources/mt19937ar.h b/CUtils/c_sources/mt19937ar.h
index 6ec8623f5f83bd70a53939449537cf11886bd2a9..16663c0f30580cbef8a040b987a38d4acb779ed9 100644
--- a/CUtils/c_sources/mt19937ar.h
+++ b/CUtils/c_sources/mt19937ar.h
@@ -41,6 +41,9 @@
    email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
 */
 
+#ifndef _MT19937AR_H
+#define _MT19937AR_H
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -68,19 +71,25 @@ unsigned long genrand_int32_mt(mt19937ar_t* priv);
 
 /* generates a random number on [0,0x7fffffff]-interval */
 long genrand_int31(void);
+long genrand_int31_mt(mt19937ar_t* priv);
 
 /* generates a random number on [0,1]-real-interval */
 double genrand_real1(void);
+double genrand_real1_mt(mt19937ar_t* priv);
 
 /* generates a random number on [0,1)-real-interval */
 double genrand_real2(void);
+double genrand_real2_mt(mt19937ar_t* priv);
 
 /* generates a random number on (0,1)-real-interval */
 double genrand_real3(void);
+double genrand_real3_mt(mt19937ar_t* priv);
 
 /* generates a random number on [0,1) with 53-bit resolution*/
 double genrand_res53(void);
+double genrand_res53_mt(mt19937ar_t* priv);
 
 #ifdef __cplusplus
 }
 #endif
+#endif