Index: trunk/SophyaLib/BaseTools/Makefile
===================================================================
--- trunk/SophyaLib/BaseTools/Makefile	(revision 3601)
+++ trunk/SophyaLib/BaseTools/Makefile	(revision 3602)
@@ -4,8 +4,8 @@
  
 clean:
-	 rm -f $(SOPHYAOBJP)datatype.o $(SOPHYAOBJP)dvlist.o $(SOPHYAOBJP)fiondblock.o $(SOPHYAOBJP)mutyv.o $(SOPHYAOBJP)ndatablock.o $(SOPHYAOBJP)pexceptions.o $(SOPHYAOBJP)ppersist.o $(SOPHYAOBJP)ppfbinstream.o $(SOPHYAOBJP)ppftpointerio.o $(SOPHYAOBJP)rawstream.o $(SOPHYAOBJP)sophyainit.o $(SOPHYAOBJP)stsrand.o $(SOPHYAOBJP)timestamp.o $(SOPHYAOBJP)gnumd5.o $(SOPHYAOBJP)srandgen.o
+	 rm -f $(SOPHYAOBJP)datatype.o $(SOPHYAOBJP)dvlist.o $(SOPHYAOBJP)fiondblock.o $(SOPHYAOBJP)mutyv.o $(SOPHYAOBJP)ndatablock.o $(SOPHYAOBJP)pexceptions.o $(SOPHYAOBJP)ppersist.o $(SOPHYAOBJP)ppfbinstream.o $(SOPHYAOBJP)ppftpointerio.o $(SOPHYAOBJP)randfmt.o $(SOPHYAOBJP)randinterf.o $(SOPHYAOBJP)randr48.o $(SOPHYAOBJP)rawstream.o $(SOPHYAOBJP)sophyainit.o $(SOPHYAOBJP)stsrand.o $(SOPHYAOBJP)timestamp.o $(SOPHYAOBJP)dSFMT.o $(SOPHYAOBJP)gnumd5.o $(SOPHYAOBJP)srandgen.o
 	rm -f $(SOPHYALIBP)libBaseTools.a
  
-$(SOPHYALIBP)libBaseTools.a : $(SOPHYAOBJP)datatype.o $(SOPHYAOBJP)dvlist.o $(SOPHYAOBJP)fiondblock.o $(SOPHYAOBJP)gnumd5.o $(SOPHYAOBJP)mutyv.o $(SOPHYAOBJP)ndatablock.o $(SOPHYAOBJP)pexceptions.o $(SOPHYAOBJP)ppersist.o $(SOPHYAOBJP)ppfbinstream.o $(SOPHYAOBJP)ppftpointerio.o $(SOPHYAOBJP)rawstream.o $(SOPHYAOBJP)sophyainit.o $(SOPHYAOBJP)srandgen.o $(SOPHYAOBJP)stsrand.o $(SOPHYAOBJP)timestamp.o
+$(SOPHYALIBP)libBaseTools.a : $(SOPHYAOBJP)datatype.o $(SOPHYAOBJP)dSFMT.o $(SOPHYAOBJP)dvlist.o $(SOPHYAOBJP)fiondblock.o $(SOPHYAOBJP)gnumd5.o $(SOPHYAOBJP)mutyv.o $(SOPHYAOBJP)ndatablock.o $(SOPHYAOBJP)pexceptions.o $(SOPHYAOBJP)ppersist.o $(SOPHYAOBJP)ppfbinstream.o $(SOPHYAOBJP)ppftpointerio.o $(SOPHYAOBJP)randfmt.o $(SOPHYAOBJP)randinterf.o $(SOPHYAOBJP)randr48.o $(SOPHYAOBJP)rawstream.o $(SOPHYAOBJP)sophyainit.o $(SOPHYAOBJP)srandgen.o $(SOPHYAOBJP)stsrand.o $(SOPHYAOBJP)timestamp.o
 	$(AR) $(ARFLAGS) $@ $? 
 	touch $(SOPHYAINCP)/SophyaConfInfo/libsophya.objlist
@@ -59,4 +59,21 @@
 	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  ppftpointerio.cc
  
+$(SOPHYAOBJP)randfmt.o: randfmt.cc randfmt.h \
+  $(SOPHYAINCP)machdefs.h objfio.h anydataobj.h \
+  ppersist.h pexceptions.h gnumd5.h ppfbinstream.h rawstream.h \
+  ppfnametag.h dsfmtflags.h dSFMT.h randinterf.h
+	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  randfmt.cc
+ 
+$(SOPHYAOBJP)randinterf.o: randinterf.cc sopnamsp.h \
+  $(SOPHYAINCP)machdefs.h pexceptions.h \
+  randinterf.h anydataobj.h
+	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  randinterf.cc
+ 
+$(SOPHYAOBJP)randr48.o: randr48.cc $(SOPHYAINCP)machdefs.h \
+  thsafeop.h pexceptions.h fiondblock.h ppersist.h gnumd5.h \
+  ppfbinstream.h rawstream.h ppfnametag.h ndatablock.h anydataobj.h \
+  randr48.h objfio.h randinterf.h
+	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  randr48.cc
+ 
 $(SOPHYAOBJP)rawstream.o: rawstream.cc sopnamsp.h rawstream.h \
   $(SOPHYAINCP)machdefs.h pexceptions.h
@@ -67,5 +84,6 @@
   pexceptions.h ppersist.h gnumd5.h ppfbinstream.h rawstream.h \
   ppfnametag.h fiondblock.h ndatablock.h anydataobj.h dvlist.h objfio.h \
-  mutyv.h timestamp.h stsrand.h fiosegdb.h segdatablock.h thsafeop.h \
+  mutyv.h timestamp.h stsrand.h randr48.h randinterf.h randfmt.h \
+  dsfmtflags.h dSFMT.h fiosegdb.h segdatablock.h thsafeop.h \
   ppftpointerio.h ppfwrapstlv.h sversion.h
 	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  sophyainit.cc
@@ -82,4 +100,7 @@
 	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  timestamp.cc
  
+$(SOPHYAOBJP)dSFMT.o: dSFMT.c dsfmtflags.h dSFMT-params.h dSFMT.h dSFMT-params19937.h
+	$(CCOMPILE) -o $@  dSFMT.c
+ 
 $(SOPHYAOBJP)gnumd5.o: gnumd5.c $(SOPHYAINCP)machdefs.h \
   bithelp.h gnumd5.h
Index: trunk/SophyaLib/BaseTools/basetools.h
===================================================================
--- trunk/SophyaLib/BaseTools/basetools.h	(revision 3601)
+++ trunk/SophyaLib/BaseTools/basetools.h	(revision 3602)
@@ -23,4 +23,6 @@
 #include "stsrand.h"
 #include "datatype.h"
+#include "randfmt.h"
+#include "randr48.h"
 
 #endif
Index: trunk/SophyaLib/BaseTools/dSFMT-params.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params.h	(revision 3602)
@@ -0,0 +1,87 @@
+#ifndef DSFMT_PARAMS_H
+#define DSFMT_PARAMS_H
+
+#include "dSFMT.h"
+
+/*----------------------
+  the parameters of DSFMT
+  following definitions are in dSFMT-paramsXXXX.h file.
+  ----------------------*/
+/** the pick up position of the array.
+#define DSFMT_POS1 122 
+*/
+
+/** the parameter of shift left as four 32-bit registers.
+#define DSFMT_SL1 18
+ */
+
+/** the parameter of shift right as four 32-bit registers.
+#define DSFMT_SR1 12
+*/
+
+/** A bitmask, used in the recursion.  These parameters are introduced
+ * to break symmetry of SIMD.
+#define DSFMT_MSK1 (uint64_t)0xdfffffefULL
+#define DSFMT_MSK2 (uint64_t)0xddfecb7fULL
+*/
+
+/** These definitions are part of a 128-bit period certification vector.
+#define DSFMT_PCV1	UINT64_C(0x00000001)
+#define DSFMT_PCV2	UINT64_C(0x00000000)
+*/
+
+#define DSFMT_LOW_MASK  UINT64_C(0x000FFFFFFFFFFFFF)
+#define DSFMT_HIGH_CONST UINT64_C(0x3FF0000000000000)
+#define DSFMT_SR	12
+
+/* for sse2 */
+#if defined(HAVE_SSE2)
+  #define SSE2_SHUFF 0x1b
+#elif defined(HAVE_ALTIVEC)
+  #if defined(__APPLE__)  /* For OSX */
+    #define ALTI_SR (vector unsigned char)(4)
+    #define ALTI_SR_PERM \
+        (vector unsigned char)(15,0,1,2,3,4,5,6,15,8,9,10,11,12,13,14)
+    #define ALTI_SR_MSK \
+        (vector unsigned int)(0x000fffffU,0xffffffffU,0x000fffffU,0xffffffffU)
+    #define ALTI_PERM \
+        (vector unsigned char)(12,13,14,15,8,9,10,11,4,5,6,7,0,1,2,3)
+  #else
+    #define ALTI_SR      {4}
+    #define ALTI_SR_PERM {15,0,1,2,3,4,5,6,15,8,9,10,11,12,13,14}
+    #define ALTI_SR_MSK  {0x000fffffU,0xffffffffU,0x000fffffU,0xffffffffU}
+    #define ALTI_PERM    {12,13,14,15,8,9,10,11,4,5,6,7,0,1,2,3}
+  #endif
+#endif
+
+#if DSFMT_MEXP == 521
+  #include "dSFMT-params521.h"
+#elif DSFMT_MEXP == 1279
+  #include "dSFMT-params1279.h"
+#elif DSFMT_MEXP == 2203
+  #include "dSFMT-params2203.h"
+#elif DSFMT_MEXP == 4253
+  #include "dSFMT-params4253.h"
+#elif DSFMT_MEXP == 11213
+  #include "dSFMT-params11213.h"
+#elif DSFMT_MEXP == 19937
+  #include "dSFMT-params19937.h"
+#elif DSFMT_MEXP == 44497
+  #include "dSFMT-params44497.h"
+#elif DSFMT_MEXP == 86243
+  #include "dSFMT-params86243.h"
+#elif DSFMT_MEXP == 132049
+  #include "dSFMT-params132049.h"
+#elif DSFMT_MEXP == 216091
+  #include "dSFMT-params216091.h"
+#else
+#ifdef __GNUC__
+  #error "DSFMT_MEXP is not valid."
+  #undef DSFMT_MEXP
+#else
+  #undef DSFMT_MEXP
+#endif
+
+#endif
+
+#endif /* DSFMT_PARAMS_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params11213.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params11213.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params11213.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS11213_H
+#define DSFMT_PARAMS11213_H
+
+/* #define DSFMT_N	107 */
+/* #define DSFMT_MAXDEGREE	11256 */
+#define DSFMT_POS1	37
+#define DSFMT_SL1	19
+#define DSFMT_MSK1	UINT64_C(0x000ffffffdf7fffd)
+#define DSFMT_MSK2	UINT64_C(0x000dfffffff6bfff)
+#define DSFMT_MSK32_1	0x000fffffU
+#define DSFMT_MSK32_2	0xfdf7fffdU
+#define DSFMT_MSK32_3	0x000dffffU
+#define DSFMT_MSK32_4	0xfff6bfffU
+#define DSFMT_FIX1	UINT64_C(0xd0ef7b7c75b06793)
+#define DSFMT_FIX2	UINT64_C(0x9c50ff4caae0a641)
+#define DSFMT_PCV1	UINT64_C(0x8234c51207c80000)
+#define DSFMT_PCV2	UINT64_C(0x0000000000000001)
+#define DSFMT_IDSTR	"dSFMT2-11213:37-19:ffffffdf7fffd-dfffffff6bfff"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(3, 3, 3, 3)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{3, 3, 3, 3}
+    #define ALTI_SL1_PERM \
+	{2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS11213_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params1279.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params1279.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params1279.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS1279_H
+#define DSFMT_PARAMS1279_H
+
+/* #define DSFMT_N	12 */
+/* #define DSFMT_MAXDEGREE	1376 */
+#define DSFMT_POS1	9
+#define DSFMT_SL1	19
+#define DSFMT_MSK1	UINT64_C(0x000efff7ffddffee)
+#define DSFMT_MSK2	UINT64_C(0x000fbffffff77fff)
+#define DSFMT_MSK32_1	0x000efff7U
+#define DSFMT_MSK32_2	0xffddffeeU
+#define DSFMT_MSK32_3	0x000fbfffU
+#define DSFMT_MSK32_4	0xfff77fffU
+#define DSFMT_FIX1	UINT64_C(0xb66627623d1a31be)
+#define DSFMT_FIX2	UINT64_C(0x04b6c51147b6109b)
+#define DSFMT_PCV1	UINT64_C(0x7049f2da382a6aeb)
+#define DSFMT_PCV2	UINT64_C(0xde4ca84a40000001)
+#define DSFMT_IDSTR	"dSFMT2-1279:9-19:efff7ffddffee-fbffffff77fff"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(3, 3, 3, 3)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{3, 3, 3, 3}
+    #define ALTI_SL1_PERM \
+	{2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS1279_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params132049.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params132049.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params132049.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS132049_H
+#define DSFMT_PARAMS132049_H
+
+/* #define DSFMT_N	1269 */
+/* #define DSFMT_MAXDEGREE	132104 */
+#define DSFMT_POS1	371
+#define DSFMT_SL1	23
+#define DSFMT_MSK1	UINT64_C(0x000fb9f4eff4bf77)
+#define DSFMT_MSK2	UINT64_C(0x000fffffbfefff37)
+#define DSFMT_MSK32_1	0x000fb9f4U
+#define DSFMT_MSK32_2	0xeff4bf77U
+#define DSFMT_MSK32_3	0x000fffffU
+#define DSFMT_MSK32_4	0xbfefff37U
+#define DSFMT_FIX1	UINT64_C(0x4ce24c0e4e234f3b)
+#define DSFMT_FIX2	UINT64_C(0x62612409b5665c2d)
+#define DSFMT_PCV1	UINT64_C(0x181232889145d000)
+#define DSFMT_PCV2	UINT64_C(0x0000000000000001)
+#define DSFMT_IDSTR	"dSFMT2-132049:371-23:fb9f4eff4bf77-fffffbfefff37"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(7, 7, 7, 7)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xff800000U,0xffffffffU,0xff800000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{7, 7, 7, 7}
+    #define ALTI_SL1_PERM \
+	{2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xff800000U,0xffffffffU,0xff800000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS132049_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params19937.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params19937.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params19937.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS19937_H
+#define DSFMT_PARAMS19937_H
+
+/* #define DSFMT_N	191 */
+/* #define DSFMT_MAXDEGREE	19992 */
+#define DSFMT_POS1	117
+#define DSFMT_SL1	19
+#define DSFMT_MSK1	UINT64_C(0x000ffafffffffb3f)
+#define DSFMT_MSK2	UINT64_C(0x000ffdfffc90fffd)
+#define DSFMT_MSK32_1	0x000ffaffU
+#define DSFMT_MSK32_2	0xfffffb3fU
+#define DSFMT_MSK32_3	0x000ffdffU
+#define DSFMT_MSK32_4	0xfc90fffdU
+#define DSFMT_FIX1	UINT64_C(0x90014964b32f4329)
+#define DSFMT_FIX2	UINT64_C(0x3b8d12ac548a7c7a)
+#define DSFMT_PCV1	UINT64_C(0x3d84e1ac0dc82880)
+#define DSFMT_PCV2	UINT64_C(0x0000000000000001)
+#define DSFMT_IDSTR	"dSFMT2-19937:117-19:ffafffffffb3f-ffdfffc90fffd"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(3, 3, 3, 3)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{3, 3, 3, 3}
+    #define ALTI_SL1_PERM \
+	{2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS19937_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params216091.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params216091.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params216091.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS216091_H
+#define DSFMT_PARAMS216091_H
+
+/* #define DSFMT_N	2077 */
+/* #define DSFMT_MAXDEGREE	216136 */
+#define DSFMT_POS1	1890
+#define DSFMT_SL1	23
+#define DSFMT_MSK1	UINT64_C(0x000bf7df7fefcfff)
+#define DSFMT_MSK2	UINT64_C(0x000e7ffffef737ff)
+#define DSFMT_MSK32_1	0x000bf7dfU
+#define DSFMT_MSK32_2	0x7fefcfffU
+#define DSFMT_MSK32_3	0x000e7fffU
+#define DSFMT_MSK32_4	0xfef737ffU
+#define DSFMT_FIX1	UINT64_C(0xd7f95a04764c27d7)
+#define DSFMT_FIX2	UINT64_C(0x6a483861810bebc2)
+#define DSFMT_PCV1	UINT64_C(0x3af0a8f3d5600000)
+#define DSFMT_PCV2	UINT64_C(0x0000000000000001)
+#define DSFMT_IDSTR	"dSFMT2-216091:1890-23:bf7df7fefcfff-e7ffffef737ff"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(7, 7, 7, 7)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xff800000U,0xffffffffU,0xff800000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{7, 7, 7, 7}
+    #define ALTI_SL1_PERM \
+	{2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xff800000U,0xffffffffU,0xff800000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS216091_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params2203.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params2203.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params2203.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS2203_H
+#define DSFMT_PARAMS2203_H
+
+/* #define DSFMT_N	20 */
+/* #define DSFMT_MAXDEGREE	2208 */
+#define DSFMT_POS1	7
+#define DSFMT_SL1	19
+#define DSFMT_MSK1	UINT64_C(0x000fdffff5edbfff)
+#define DSFMT_MSK2	UINT64_C(0x000f77fffffffbfe)
+#define DSFMT_MSK32_1	0x000fdfffU
+#define DSFMT_MSK32_2	0xf5edbfffU
+#define DSFMT_MSK32_3	0x000f77ffU
+#define DSFMT_MSK32_4	0xfffffbfeU
+#define DSFMT_FIX1	UINT64_C(0xb14e907a39338485)
+#define DSFMT_FIX2	UINT64_C(0xf98f0735c637ef90)
+#define DSFMT_PCV1	UINT64_C(0x8000000000000000)
+#define DSFMT_PCV2	UINT64_C(0x0000000000000001)
+#define DSFMT_IDSTR	"dSFMT2-2203:7-19:fdffff5edbfff-f77fffffffbfe"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(3, 3, 3, 3)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{3, 3, 3, 3}
+    #define ALTI_SL1_PERM \
+	{2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS2203_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params4253.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params4253.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params4253.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS4253_H
+#define DSFMT_PARAMS4253_H
+
+/* #define DSFMT_N	40 */
+/* #define DSFMT_MAXDEGREE	4288 */
+#define DSFMT_POS1	19
+#define DSFMT_SL1	19
+#define DSFMT_MSK1	UINT64_C(0x0007b7fffef5feff)
+#define DSFMT_MSK2	UINT64_C(0x000ffdffeffefbfc)
+#define DSFMT_MSK32_1	0x0007b7ffU
+#define DSFMT_MSK32_2	0xfef5feffU
+#define DSFMT_MSK32_3	0x000ffdffU
+#define DSFMT_MSK32_4	0xeffefbfcU
+#define DSFMT_FIX1	UINT64_C(0x80901b5fd7a11c65)
+#define DSFMT_FIX2	UINT64_C(0x5a63ff0e7cb0ba74)
+#define DSFMT_PCV1	UINT64_C(0x1ad277be12000000)
+#define DSFMT_PCV2	UINT64_C(0x0000000000000001)
+#define DSFMT_IDSTR	"dSFMT2-4253:19-19:7b7fffef5feff-ffdffeffefbfc"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(3, 3, 3, 3)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{3, 3, 3, 3}
+    #define ALTI_SL1_PERM \
+	{2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS4253_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params44497.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params44497.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params44497.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS44497_H
+#define DSFMT_PARAMS44497_H
+
+/* #define DSFMT_N	427 */
+/* #define DSFMT_MAXDEGREE	44536 */
+#define DSFMT_POS1	304
+#define DSFMT_SL1	19
+#define DSFMT_MSK1	UINT64_C(0x000ff6dfffffffef)
+#define DSFMT_MSK2	UINT64_C(0x0007ffdddeefff6f)
+#define DSFMT_MSK32_1	0x000ff6dfU
+#define DSFMT_MSK32_2	0xffffffefU
+#define DSFMT_MSK32_3	0x0007ffddU
+#define DSFMT_MSK32_4	0xdeefff6fU
+#define DSFMT_FIX1	UINT64_C(0x75d910f235f6e10e)
+#define DSFMT_FIX2	UINT64_C(0x7b32158aedc8e969)
+#define DSFMT_PCV1	UINT64_C(0x4c3356b2a0000000)
+#define DSFMT_PCV2	UINT64_C(0x0000000000000001)
+#define DSFMT_IDSTR	"dSFMT2-44497:304-19:ff6dfffffffef-7ffdddeefff6f"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(3, 3, 3, 3)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{3, 3, 3, 3}
+    #define ALTI_SL1_PERM \
+	{2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS44497_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params521.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params521.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params521.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS521_H
+#define DSFMT_PARAMS521_H
+
+/* #define DSFMT_N	4 */
+/* #define DSFMT_MAXDEGREE	544 */
+#define DSFMT_POS1	3
+#define DSFMT_SL1	25
+#define DSFMT_MSK1	UINT64_C(0x000fbfefff77efff)
+#define DSFMT_MSK2	UINT64_C(0x000ffeebfbdfbfdf)
+#define DSFMT_MSK32_1	0x000fbfefU
+#define DSFMT_MSK32_2	0xff77efffU
+#define DSFMT_MSK32_3	0x000ffeebU
+#define DSFMT_MSK32_4	0xfbdfbfdfU
+#define DSFMT_FIX1	UINT64_C(0xcfb393d661638469)
+#define DSFMT_FIX2	UINT64_C(0xc166867883ae2adb)
+#define DSFMT_PCV1	UINT64_C(0xccaa588000000000)
+#define DSFMT_PCV2	UINT64_C(0x0000000000000001)
+#define DSFMT_IDSTR	"dSFMT2-521:3-25:fbfefff77efff-ffeebfbdfbfdf"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(1, 1, 1, 1)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xfe000000U,0xffffffffU,0xfe000000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{1, 1, 1, 1}
+    #define ALTI_SL1_PERM \
+	{3,4,5,6,7,29,29,29,11,12,13,14,15,0,1,2}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xfe000000U,0xffffffffU,0xfe000000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS521_H */
Index: trunk/SophyaLib/BaseTools/dSFMT-params86243.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT-params86243.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT-params86243.h	(revision 3602)
@@ -0,0 +1,40 @@
+#ifndef DSFMT_PARAMS86243_H
+#define DSFMT_PARAMS86243_H
+
+/* #define DSFMT_N	829 */
+/* #define DSFMT_MAXDEGREE	86344 */
+#define DSFMT_POS1	231
+#define DSFMT_SL1	13
+#define DSFMT_MSK1	UINT64_C(0x000ffedff6ffffdf)
+#define DSFMT_MSK2	UINT64_C(0x000ffff7fdffff7e)
+#define DSFMT_MSK32_1	0x000ffedfU
+#define DSFMT_MSK32_2	0xf6ffffdfU
+#define DSFMT_MSK32_3	0x000ffff7U
+#define DSFMT_MSK32_4	0xfdffff7eU
+#define DSFMT_FIX1	UINT64_C(0x1d553e776b975e68)
+#define DSFMT_FIX2	UINT64_C(0x648faadf1416bf91)
+#define DSFMT_PCV1	UINT64_C(0x5f2cd03e2758a373)
+#define DSFMT_PCV2	UINT64_C(0xc0b7eb8410000001)
+#define DSFMT_IDSTR	"dSFMT2-86243:231-13:ffedff6ffffdf-ffff7fdffff7e"
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__)	/* For OSX */
+    #define ALTI_SL1 	(vector unsigned int)(5, 5, 5, 5)
+    #define ALTI_SL1_PERM \
+	(vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
+    #define ALTI_SL1_MSK \
+	(vector unsigned int)(0xffffffffU,0xffffe000U,0xffffffffU,0xffffe000U)
+    #define ALTI_MSK	(vector unsigned int)(DSFMT_MSK32_1, \
+			DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4)
+#else	/* For OTHER OSs(Linux?) */
+    #define ALTI_SL1 	{5, 5, 5, 5}
+    #define ALTI_SL1_PERM \
+	{1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
+    #define ALTI_SL1_MSK \
+	{0xffffffffU,0xffffe000U,0xffffffffU,0xffffe000U}
+    #define ALTI_MSK \
+	{DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4}
+#endif
+
+#endif /* DSFMT_PARAMS86243_H */
Index: trunk/SophyaLib/BaseTools/dSFMT.c
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT.c	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT.c	(revision 3602)
@@ -0,0 +1,735 @@
+/** 
+ * @file dSFMT.c 
+ * @brief double precision SIMD-oriented Fast Mersenne Twister (dSFMT)
+ * based on IEEE 754 format.
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * Copyright (C) 2007,2008 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University. All rights reserved.
+ *
+ * The new BSD License is applied to this software, see LICENSE.txt
+ */
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include "dsfmtflags.h"
+#include "dSFMT-params.h"
+
+/** dsfmt internal state vector */
+dsfmt_t dsfmt_global_data;
+/** dsfmt mexp for check */
+static const int dsfmt_mexp = DSFMT_MEXP;
+
+/*----------------
+  STATIC FUNCTIONS
+  ----------------*/
+inline static uint32_t ini_func1(uint32_t x);
+inline static uint32_t ini_func2(uint32_t x);
+inline static void gen_rand_array_c1o2(dsfmt_t *dsfmt, w128_t *array,
+				       int size);
+inline static void gen_rand_array_c0o1(dsfmt_t *dsfmt, w128_t *array,
+				       int size);
+inline static void gen_rand_array_o0c1(dsfmt_t *dsfmt, w128_t *array,
+				       int size);
+inline static void gen_rand_array_o0o1(dsfmt_t *dsfmt, w128_t *array,
+				       int size);
+inline static int idxof(int i);
+static void initial_mask(dsfmt_t *dsfmt);
+static void period_certification(dsfmt_t *dsfmt);
+
+#if defined(HAVE_SSE2)
+#  include <emmintrin.h>
+/** mask data for sse2 */
+static __m128i sse2_param_mask;
+/** 1 in 64bit for sse2 */
+static __m128i sse2_int_one;
+/** 2.0 double for sse2 */
+static __m128d sse2_double_two;
+/** -1.0 double for sse2 */
+static __m128d sse2_double_m_one;
+
+static void setup_const(void);
+#endif
+
+/**
+ * This function simulate a 32-bit array index overlapped to 64-bit
+ * array of LITTLE ENDIAN in BIG ENDIAN machine.
+ */
+#if defined(DSFMT_BIG_ENDIAN)
+inline static int idxof(int i) {
+    return i ^ 1;
+}
+#else
+inline static int idxof(int i) {
+    return i;
+}
+#endif
+
+/**
+ * This function represents the recursion formula.
+ * @param r output
+ * @param a a 128-bit part of the internal state array
+ * @param b a 128-bit part of the internal state array
+ * @param lung a 128-bit part of the internal state array
+ */
+#if defined(HAVE_ALTIVEC)
+inline static void do_recursion(w128_t *r, w128_t *a, w128_t * b,
+				w128_t *lung) {
+    const vector unsigned char sl1 = ALTI_SL1;
+    const vector unsigned char sl1_perm = ALTI_SL1_PERM;
+    const vector unsigned int sl1_msk = ALTI_SL1_MSK;
+    const vector unsigned char sr1 = ALTI_SR;
+    const vector unsigned char sr1_perm = ALTI_SR_PERM;
+    const vector unsigned int sr1_msk = ALTI_SR_MSK;
+    const vector unsigned char perm = ALTI_PERM;
+    const vector unsigned int msk1 = ALTI_MSK;
+    vector unsigned int w, x, y, z;
+
+    z = a->s;
+    w = lung->s;
+    x = vec_perm(w, (vector unsigned int)perm, perm);
+    y = vec_perm(z, sl1_perm, sl1_perm);
+    y = vec_sll(y, sl1);
+    y = vec_and(y, sl1_msk);
+    w = vec_xor(x, b->s);
+    w = vec_xor(w, y);
+    x = vec_perm(w, (vector unsigned int)sr1_perm, sr1_perm);
+    x = vec_srl(x, sr1);
+    x = vec_and(x, sr1_msk);
+    y = vec_and(w, msk1);
+    z = vec_xor(z, y);
+    r->s = vec_xor(z, x);
+    lung->s = w;
+}
+#elif defined(HAVE_SSE2)
+/**
+ * This function setup some constant variables for SSE2.
+ */
+static void setup_const(void) {
+    static int first = 1;
+    if (!first) {
+	return;
+    }
+    sse2_param_mask = _mm_set_epi32(DSFMT_MSK32_3, DSFMT_MSK32_4,
+				    DSFMT_MSK32_1, DSFMT_MSK32_2);
+    sse2_int_one = _mm_set_epi32(0, 1, 0, 1);
+    sse2_double_two = _mm_set_pd(2.0, 2.0);
+    sse2_double_m_one = _mm_set_pd(-1.0, -1.0);
+    first = 0;
+}
+
+/**
+ * This function represents the recursion formula.
+ * @param r output 128-bit
+ * @param a a 128-bit part of the internal state array
+ * @param b a 128-bit part of the internal state array
+ * @param d a 128-bit part of the internal state array (I/O)
+ */
+inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *u) {
+    __m128i v, w, x, y, z;
+    
+    x = a->si;
+    z = _mm_slli_epi64(x, DSFMT_SL1);
+    y = _mm_shuffle_epi32(u->si, SSE2_SHUFF);
+    z = _mm_xor_si128(z, b->si);
+    y = _mm_xor_si128(y, z);
+
+    v = _mm_srli_epi64(y, DSFMT_SR);
+    w = _mm_and_si128(y, sse2_param_mask);
+    v = _mm_xor_si128(v, x);
+    v = _mm_xor_si128(v, w);
+    r->si = v;
+    u->si = y;
+}
+#else /* standard C */
+/**
+ * This function represents the recursion formula.
+ * @param r output 128-bit
+ * @param a a 128-bit part of the internal state array
+ * @param b a 128-bit part of the internal state array
+ * @param lung a 128-bit part of the internal state array (I/O)
+ */
+inline static void do_recursion(w128_t *r, w128_t *a, w128_t * b,
+				w128_t *lung) {
+    uint64_t t0, t1, L0, L1;
+
+    t0 = a->u[0];
+    t1 = a->u[1];
+    L0 = lung->u[0];
+    L1 = lung->u[1];
+    lung->u[0] = (t0 << DSFMT_SL1) ^ (L1 >> 32) ^ (L1 << 32) ^ b->u[0];
+    lung->u[1] = (t1 << DSFMT_SL1) ^ (L0 >> 32) ^ (L0 << 32) ^ b->u[1];
+    r->u[0] = (lung->u[0] >> DSFMT_SR) ^ (lung->u[0] & DSFMT_MSK1) ^ t0;
+    r->u[1] = (lung->u[1] >> DSFMT_SR) ^ (lung->u[1] & DSFMT_MSK2) ^ t1;
+}
+#endif
+
+#if defined(HAVE_SSE2)
+/**
+ * This function converts the double precision floating point numbers which
+ * distribute uniformly in the range [1, 2) to those which distribute uniformly
+ * in the range [0, 1).
+ * @param w 128bit stracture of double precision floating point numbers (I/O)
+ */
+inline static void convert_c0o1(w128_t *w) {
+    w->sd = _mm_add_pd(w->sd, sse2_double_m_one);
+}
+
+/**
+ * This function converts the double precision floating point numbers which
+ * distribute uniformly in the range [1, 2) to those which distribute uniformly
+ * in the range (0, 1].
+ * @param w 128bit stracture of double precision floating point numbers (I/O)
+ */
+inline static void convert_o0c1(w128_t *w) {
+    w->sd = _mm_sub_pd(sse2_double_two, w->sd);
+}
+
+/**
+ * This function converts the double precision floating point numbers which
+ * distribute uniformly in the range [1, 2) to those which distribute uniformly
+ * in the range (0, 1).
+ * @param w 128bit stracture of double precision floating point numbers (I/O)
+ */
+inline static void convert_o0o1(w128_t *w) {
+    w->si = _mm_or_si128(w->si, sse2_int_one);
+    w->sd = _mm_add_pd(w->sd, sse2_double_m_one);
+}
+#else /* standard C and altivec */
+/**
+ * This function converts the double precision floating point numbers which
+ * distribute uniformly in the range [1, 2) to those which distribute uniformly
+ * in the range [0, 1).
+ * @param w 128bit stracture of double precision floating point numbers (I/O)
+ */
+inline static void convert_c0o1(w128_t *w) {
+    w->d[0] -= 1.0;
+    w->d[1] -= 1.0;
+}
+
+/**
+ * This function converts the double precision floating point numbers which
+ * distribute uniformly in the range [1, 2) to those which distribute uniformly
+ * in the range (0, 1].
+ * @param w 128bit stracture of double precision floating point numbers (I/O)
+ */
+inline static void convert_o0c1(w128_t *w) {
+    w->d[0] = 2.0 - w->d[0];
+    w->d[1] = 2.0 - w->d[1];
+}
+
+/**
+ * This function converts the double precision floating point numbers which
+ * distribute uniformly in the range [1, 2) to those which distribute uniformly
+ * in the range (0, 1).
+ * @param w 128bit stracture of double precision floating point numbers (I/O)
+ */
+inline static void convert_o0o1(w128_t *w) {
+    w->u[0] |= 1;
+    w->u[1] |= 1;
+    w->d[0] -= 1.0;
+    w->d[1] -= 1.0;
+}
+#endif
+
+/**
+ * This function fills the user-specified array with double precision
+ * floating point pseudorandom numbers of the IEEE 754 format.
+ * @param dsfmt dsfmt state vector.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.  
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+inline static void gen_rand_array_c1o2(dsfmt_t *dsfmt, w128_t *array,
+				       int size) {
+    int i, j;
+    w128_t lung;
+
+    lung = dsfmt->status[DSFMT_N];
+    do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1],
+		 &lung);
+    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
+	do_recursion(&array[i], &dsfmt->status[i],
+		     &dsfmt->status[i + DSFMT_POS1], &lung);
+    }
+    for (; i < DSFMT_N; i++) {
+	do_recursion(&array[i], &dsfmt->status[i],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+    }
+    for (; i < size - DSFMT_N; i++) {
+	do_recursion(&array[i], &array[i - DSFMT_N],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+    }
+    for (j = 0; j < 2 * DSFMT_N - size; j++) {
+	dsfmt->status[j] = array[j + size - DSFMT_N];
+    }
+    for (; i < size; i++, j++) {
+	do_recursion(&array[i], &array[i - DSFMT_N],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+	dsfmt->status[j] = array[i];
+    }
+    dsfmt->status[DSFMT_N] = lung;
+}
+
+/**
+ * This function fills the user-specified array with double precision
+ * floating point pseudorandom numbers of the IEEE 754 format.
+ * @param dsfmt dsfmt state vector.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.  
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+inline static void gen_rand_array_c0o1(dsfmt_t *dsfmt, w128_t *array,
+				       int size) {
+    int i, j;
+    w128_t lung;
+
+    lung = dsfmt->status[DSFMT_N];
+    do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1],
+		 &lung);
+    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
+	do_recursion(&array[i], &dsfmt->status[i],
+		     &dsfmt->status[i + DSFMT_POS1], &lung);
+    }
+    for (; i < DSFMT_N; i++) {
+	do_recursion(&array[i], &dsfmt->status[i],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+    }
+    for (; i < size - DSFMT_N; i++) {
+	do_recursion(&array[i], &array[i - DSFMT_N],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+	convert_c0o1(&array[i - DSFMT_N]);
+    }
+    for (j = 0; j < 2 * DSFMT_N - size; j++) {
+	dsfmt->status[j] = array[j + size - DSFMT_N];
+    }
+    for (; i < size; i++, j++) {
+	do_recursion(&array[i], &array[i - DSFMT_N],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+	dsfmt->status[j] = array[i];
+	convert_c0o1(&array[i - DSFMT_N]);
+    }
+    for (i = size - DSFMT_N; i < size; i++) {
+	convert_c0o1(&array[i]);
+    }
+    dsfmt->status[DSFMT_N] = lung;
+}
+
+/**
+ * This function fills the user-specified array with double precision
+ * floating point pseudorandom numbers of the IEEE 754 format.
+ * @param dsfmt dsfmt state vector.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.  
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+inline static void gen_rand_array_o0o1(dsfmt_t *dsfmt, w128_t *array,
+				       int size) {
+    int i, j;
+    w128_t lung;
+
+    lung = dsfmt->status[DSFMT_N];
+    do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1],
+		 &lung);
+    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
+	do_recursion(&array[i], &dsfmt->status[i],
+		     &dsfmt->status[i + DSFMT_POS1], &lung);
+    }
+    for (; i < DSFMT_N; i++) {
+	do_recursion(&array[i], &dsfmt->status[i],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+    }
+    for (; i < size - DSFMT_N; i++) {
+	do_recursion(&array[i], &array[i - DSFMT_N],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+	convert_o0o1(&array[i - DSFMT_N]);
+    }
+    for (j = 0; j < 2 * DSFMT_N - size; j++) {
+	dsfmt->status[j] = array[j + size - DSFMT_N];
+    }
+    for (; i < size; i++, j++) {
+	do_recursion(&array[i], &array[i - DSFMT_N],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+	dsfmt->status[j] = array[i];
+	convert_o0o1(&array[i - DSFMT_N]);
+    }
+    for (i = size - DSFMT_N; i < size; i++) {
+	convert_o0o1(&array[i]);
+    }
+    dsfmt->status[DSFMT_N] = lung;
+}
+
+/**
+ * This function fills the user-specified array with double precision
+ * floating point pseudorandom numbers of the IEEE 754 format.
+ * @param dsfmt dsfmt state vector.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.  
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+inline static void gen_rand_array_o0c1(dsfmt_t *dsfmt, w128_t *array,
+				       int size) {
+    int i, j;
+    w128_t lung;
+
+    lung = dsfmt->status[DSFMT_N];
+    do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1],
+		 &lung);
+    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
+	do_recursion(&array[i], &dsfmt->status[i],
+		     &dsfmt->status[i + DSFMT_POS1], &lung);
+    }
+    for (; i < DSFMT_N; i++) {
+	do_recursion(&array[i], &dsfmt->status[i],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+    }
+    for (; i < size - DSFMT_N; i++) {
+	do_recursion(&array[i], &array[i - DSFMT_N],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+	convert_o0c1(&array[i - DSFMT_N]);
+    }
+    for (j = 0; j < 2 * DSFMT_N - size; j++) {
+	dsfmt->status[j] = array[j + size - DSFMT_N];
+    }
+    for (; i < size; i++, j++) {
+	do_recursion(&array[i], &array[i - DSFMT_N],
+		     &array[i + DSFMT_POS1 - DSFMT_N], &lung);
+	dsfmt->status[j] = array[i];
+	convert_o0c1(&array[i - DSFMT_N]);
+    }
+    for (i = size - DSFMT_N; i < size; i++) {
+	convert_o0c1(&array[i]);
+    }
+    dsfmt->status[DSFMT_N] = lung;
+}
+
+/**
+ * This function represents a function used in the initialization
+ * by init_by_array
+ * @param x 32-bit integer
+ * @return 32-bit integer
+ */
+static uint32_t ini_func1(uint32_t x) {
+    return (x ^ (x >> 27)) * (uint32_t)1664525UL;
+}
+
+/**
+ * This function represents a function used in the initialization
+ * by init_by_array
+ * @param x 32-bit integer
+ * @return 32-bit integer
+ */
+static uint32_t ini_func2(uint32_t x) {
+    return (x ^ (x >> 27)) * (uint32_t)1566083941UL;
+}
+
+/**
+ * This function initializes the internal state array to fit the IEEE
+ * 754 format.
+ * @param dsfmt dsfmt state vector.
+ */
+static void initial_mask(dsfmt_t *dsfmt) {
+    int i;
+    uint64_t *psfmt;
+
+    psfmt = &dsfmt->status[0].u[0];
+    for (i = 0; i < DSFMT_N * 2; i++) {
+        psfmt[i] = (psfmt[i] & DSFMT_LOW_MASK) | DSFMT_HIGH_CONST;
+    }
+}
+
+/**
+ * This function certificate the period of 2^{SFMT_MEXP}-1.
+ * @param dsfmt dsfmt state vector.
+ */
+static void period_certification(dsfmt_t *dsfmt) {
+    uint64_t pcv[2] = {DSFMT_PCV1, DSFMT_PCV2};
+    uint64_t tmp[2];
+    uint64_t inner;
+    int i;
+#if (DSFMT_PCV2 & 1) != 1
+    int j;
+    uint64_t work;
+#endif
+
+    tmp[0] = (dsfmt->status[DSFMT_N].u[0] ^ DSFMT_FIX1);
+    tmp[1] = (dsfmt->status[DSFMT_N].u[1] ^ DSFMT_FIX2);
+
+    inner = tmp[0] & pcv[0];
+    inner ^= tmp[1] & pcv[1];
+    for (i = 32; i > 0; i >>= 1) {
+        inner ^= inner >> i;
+    }
+    inner &= 1;
+    /* check OK */
+    if (inner == 1) {
+	return;
+    }
+    /* check NG, and modification */
+#if (DSFMT_PCV2 & 1) == 1
+    dsfmt->status[DSFMT_N].u[1] ^= 1;
+#else
+    for (i = 1; i >= 0; i--) {
+	work = 1;
+	for (j = 0; j < 64; j++) {
+	    if ((work & pcv[i]) != 0) {
+		dsfmt->status[DSFMT_N].u[i] ^= work;
+		return;
+	    }
+	    work = work << 1;
+	}
+    }
+#endif
+    return;
+}
+
+/*----------------
+  PUBLIC FUNCTIONS
+  ----------------*/
+/**
+ * This function returns the identification string.  The string shows
+ * the Mersenne exponent, and all parameters of this generator.
+ * @return id string.
+ */
+const char *dsfmt_get_idstring(void) {
+    return DSFMT_IDSTR;
+}
+
+/**
+ * This function returns the minimum size of array used for \b
+ * fill_array functions.
+ * @return minimum size of array used for fill_array functions.
+ */
+int dsfmt_get_min_array_size(void) {
+    return DSFMT_N64;
+}
+
+/**
+ * This function fills the internal state array with double precision
+ * floating point pseudorandom numbers of the IEEE 754 format.
+ * @param dsfmt dsfmt state vector.
+ */
+void dsfmt_gen_rand_all(dsfmt_t *dsfmt) {
+    int i;
+    w128_t lung;
+
+    lung = dsfmt->status[DSFMT_N];
+    do_recursion(&dsfmt->status[0], &dsfmt->status[0],
+		 &dsfmt->status[DSFMT_POS1], &lung);
+    for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) {
+	do_recursion(&dsfmt->status[i], &dsfmt->status[i],
+		     &dsfmt->status[i + DSFMT_POS1], &lung);
+    }
+    for (; i < DSFMT_N; i++) {
+	do_recursion(&dsfmt->status[i], &dsfmt->status[i],
+		     &dsfmt->status[i + DSFMT_POS1 - DSFMT_N], &lung);
+    }
+    dsfmt->status[DSFMT_N] = lung;
+}
+
+/**
+ * This function generates double precision floating point
+ * pseudorandom numbers which distribute in the range [1, 2) to the
+ * specified array[] by one call. The number of pseudorandom numbers
+ * is specified by the argument \b size, which must be at least (SFMT_MEXP
+ * / 128) * 2 and a multiple of two.  The function
+ * get_min_array_size() returns this minimum size.  The generation by
+ * this function is much faster than the following fill_array_xxx functions.
+ *
+ * For initialization, init_gen_rand() or init_by_array() must be called
+ * before the first call of this function. This function can not be
+ * used after calling genrand_xxx functions, without initialization.
+ *
+ * @param dsfmt dsfmt state vector.
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.  The pointer to the array must be "aligned"
+ * (namely, must be a multiple of 16) in the SIMD version, since it
+ * refers to the address of a 128-bit integer.  In the standard C
+ * version, the pointer is arbitrary.
+ *
+ * @param size the number of 64-bit pseudorandom integers to be
+ * generated.  size must be a multiple of 2, and greater than or equal
+ * to (SFMT_MEXP / 128) * 2.
+ *
+ * @note \b memalign or \b posix_memalign is available to get aligned
+ * memory. Mac OSX doesn't have these functions, but \b malloc of OSX
+ * returns the pointer to the aligned memory block.
+ */
+void dsfmt_fill_array_close1_open2(dsfmt_t *dsfmt, double array[], int size) {
+    assert(size % 2 == 0);
+    assert(size >= DSFMT_N64);
+    gen_rand_array_c1o2(dsfmt, (w128_t *)array, size / 2);
+}
+
+/**
+ * This function generates double precision floating point
+ * pseudorandom numbers which distribute in the range (0, 1] to the
+ * specified array[] by one call. This function is the same as
+ * fill_array_close1_open2() except the distribution range.
+ *
+ * @param dsfmt dsfmt state vector.
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa fill_array_close1_open2()
+ */
+void dsfmt_fill_array_open_close(dsfmt_t *dsfmt, double array[], int size) {
+    assert(size % 2 == 0);
+    assert(size >= DSFMT_N64);
+    gen_rand_array_o0c1(dsfmt, (w128_t *)array, size / 2);
+}
+
+/**
+ * This function generates double precision floating point
+ * pseudorandom numbers which distribute in the range [0, 1) to the
+ * specified array[] by one call. This function is the same as
+ * fill_array_close1_open2() except the distribution range.
+ *
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param dsfmt dsfmt state vector.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa fill_array_close1_open2()
+ */
+void dsfmt_fill_array_close_open(dsfmt_t *dsfmt, double array[], int size) {
+    assert(size % 2 == 0);
+    assert(size >= DSFMT_N64);
+    gen_rand_array_c0o1(dsfmt, (w128_t *)array, size / 2);
+}
+
+/**
+ * This function generates double precision floating point
+ * pseudorandom numbers which distribute in the range (0, 1) to the
+ * specified array[] by one call. This function is the same as
+ * fill_array_close1_open2() except the distribution range.
+ *
+ * @param dsfmt dsfmt state vector.
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa fill_array_close1_open2()
+ */
+void dsfmt_fill_array_open_open(dsfmt_t *dsfmt, double array[], int size) {
+    assert(size % 2 == 0);
+    assert(size >= DSFMT_N64);
+    gen_rand_array_o0o1(dsfmt, (w128_t *)array, size / 2);
+}
+
+#if defined(__INTEL_COMPILER)
+#  pragma warning(disable:981)
+#endif
+/**
+ * This function initializes the internal state array with a 32-bit
+ * integer seed.
+ * @param dsfmt dsfmt state vector.
+ * @param seed a 32-bit integer used as the seed.
+ * @param mexp caller's mersenne expornent
+ */
+void dsfmt_chk_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed, int mexp) {
+    int i;
+    uint32_t *psfmt;
+
+    /* make sure caller program is compiled with the same MEXP */
+    if (mexp != dsfmt_mexp) {
+	fprintf(stderr, "DSFMT_MEXP doesn't match with dSFMT.c\n");
+	exit(1);
+    }
+    psfmt = &dsfmt->status[0].u32[0];
+    psfmt[idxof(0)] = seed;
+    for (i = 1; i < (DSFMT_N + 1) * 4; i++) {
+        psfmt[idxof(i)] = 1812433253UL 
+	    * (psfmt[idxof(i - 1)] ^ (psfmt[idxof(i - 1)] >> 30)) + i;
+    }
+    initial_mask(dsfmt);
+    period_certification(dsfmt);
+    dsfmt->idx = DSFMT_N64;
+#if defined(HAVE_SSE2)
+    setup_const();
+#endif
+}
+
+/**
+ * This function initializes the internal state array,
+ * with an array of 32-bit integers used as the seeds
+ * @param dsfmt dsfmt state vector.
+ * @param init_key the array of 32-bit integers, used as a seed.
+ * @param key_length the length of init_key.
+ * @param mexp caller's mersenne expornent
+ */
+void dsfmt_chk_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[],
+			     int key_length, int mexp) {
+    int i, j, count;
+    uint32_t r;
+    uint32_t *psfmt32;
+    int lag;
+    int mid;
+    int size = (DSFMT_N + 1) * 4;	/* pulmonary */
+
+    /* make sure caller program is compiled with the same MEXP */
+    if (mexp != dsfmt_mexp) {
+	fprintf(stderr, "DSFMT_MEXP doesn't match with dSFMT.c\n");
+	exit(1);
+    }
+    if (size >= 623) {
+	lag = 11;
+    } else if (size >= 68) {
+	lag = 7;
+    } else if (size >= 39) {
+	lag = 5;
+    } else {
+	lag = 3;
+    }
+    mid = (size - lag) / 2;
+
+    psfmt32 = &dsfmt->status[0].u32[0];
+    memset(dsfmt->status, 0x8b, sizeof(dsfmt->status));
+    if (key_length + 1 > size) {
+	count = key_length + 1;
+    } else {
+	count = size;
+    }
+    r = ini_func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid % size)] 
+		  ^ psfmt32[idxof((size - 1) % size)]);
+    psfmt32[idxof(mid % size)] += r;
+    r += key_length;
+    psfmt32[idxof((mid + lag) % size)] += r;
+    psfmt32[idxof(0)] = r;
+    count--;
+    for (i = 1, j = 0; (j < count) && (j < key_length); j++) {
+	r = ini_func1(psfmt32[idxof(i)] 
+		      ^ psfmt32[idxof((i + mid) % size)] 
+		      ^ psfmt32[idxof((i + size - 1) % size)]);
+	psfmt32[idxof((i + mid) % size)] += r;
+	r += init_key[j] + i;
+	psfmt32[idxof((i + mid + lag) % size)] += r;
+	psfmt32[idxof(i)] = r;
+	i = (i + 1) % size;
+    }
+    for (; j < count; j++) {
+	r = ini_func1(psfmt32[idxof(i)] 
+		      ^ psfmt32[idxof((i + mid) % size)] 
+		      ^ psfmt32[idxof((i + size - 1) % size)]);
+	psfmt32[idxof((i + mid) % size)] += r;
+	r += i;
+	psfmt32[idxof((i + mid + lag) % size)] += r;
+	psfmt32[idxof(i)] = r;
+	i = (i + 1) % size;
+    }
+    for (j = 0; j < size; j++) {
+	r = ini_func2(psfmt32[idxof(i)] 
+		      + psfmt32[idxof((i + mid) % size)] 
+		      + psfmt32[idxof((i + size - 1) % size)]);
+	psfmt32[idxof((i + mid) % size)] ^= r;
+	r -= i;
+	psfmt32[idxof((i + mid + lag) % size)] ^= r;
+	psfmt32[idxof(i)] = r;
+	i = (i + 1) % size;
+    }
+    initial_mask(dsfmt);
+    period_certification(dsfmt);
+    dsfmt->idx = DSFMT_N64;
+#if defined(HAVE_SSE2)
+    setup_const();
+#endif
+}
+#if defined(__INTEL_COMPILER)
+#  pragma warning(default:981)
+#endif
Index: trunk/SophyaLib/BaseTools/dSFMT.h
===================================================================
--- trunk/SophyaLib/BaseTools/dSFMT.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dSFMT.h	(revision 3602)
@@ -0,0 +1,632 @@
+/** 
+ * @file dSFMT.h 
+ *
+ * @brief double precision SIMD oriented Fast Mersenne Twister(dSFMT)
+ * pseudorandom number generator based on IEEE 754 format.
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * Copyright (C) 2007, 2008 Mutsuo Saito, Makoto Matsumoto and
+ * Hiroshima University. All rights reserved.
+ *
+ * The new BSD License is applied to this software.
+ * see LICENSE.txt
+ *
+ * @note We assume that your system has inttypes.h.  If your system
+ * doesn't have inttypes.h, you have to typedef uint32_t and uint64_t,
+ * and you have to define PRIu64 and PRIx64 in this file as follows:
+ * @verbatim
+ typedef unsigned int uint32_t
+ typedef unsigned long long uint64_t  
+ #define PRIu64 "llu"
+ #define PRIx64 "llx"
+@endverbatim
+ * uint32_t must be exactly 32-bit unsigned integer type (no more, no
+ * less), and uint64_t must be exactly 64-bit unsigned integer type.
+ * PRIu64 and PRIx64 are used for printf function to print 64-bit
+ * unsigned int and 64-bit unsigned int in hexadecimal format.
+ */
+
+#ifndef DSFMT_H
+#define DSFMT_H
+
+#include <stdio.h>
+#include <assert.h>
+
+#if !defined(DSFMT_MEXP)
+#ifdef __GNUC__
+  #warning "DSFMT_MEXP is not defined. I assume DSFMT_MEXP is 19937."
+#endif
+  #define DSFMT_MEXP 19937
+#endif
+/*-----------------
+  BASIC DEFINITIONS
+  -----------------*/
+/* Mersenne Exponent. The period of the sequence 
+ *  is a multiple of 2^DSFMT_MEXP-1.
+ * #define DSFMT_MEXP 19937 */
+/** DSFMT generator has an internal state array of 128-bit integers,
+ * and N is its size. */
+#define DSFMT_N ((DSFMT_MEXP - 128) / 104 + 1)
+/** N32 is the size of internal state array when regarded as an array
+ * of 32-bit integers.*/
+#define DSFMT_N32 (DSFMT_N * 4)
+/** N64 is the size of internal state array when regarded as an array
+ * of 64-bit integers.*/
+#define DSFMT_N64 (DSFMT_N * 2)
+
+#if !defined(DSFMT_BIG_ENDIAN)
+#  if defined(__BYTE_ORDER) && defined(__BIG_ENDIAN)
+#    if __BYTE_ORDER == __BIG_ENDIAN
+#      define DSFMT_BIG_ENDIAN 1
+#    endif
+#  elif defined(_BYTE_ORDER) && defined(_BIG_ENDIAN)
+#    if _BYTE_ORDER == _BIG_ENDIAN
+#      define DSFMT_BIG_ENDIAN 1
+#    endif
+#  elif defined(__BYTE_ORDER__) && defined(__BIG_ENDIAN__)
+#    if __BYTE_ORDER__ == __BIG_ENDIAN__
+#      define DSFMT_BIG_ENDIAN 1
+#    endif
+#  elif defined(BYTE_ORDER) && defined(BIG_ENDIAN)
+#    if BYTE_ORDER == BIG_ENDIAN
+#      define DSFMT_BIG_ENDIAN 1
+#    endif
+#  elif defined(__BIG_ENDIAN) || defined(_BIG_ENDIAN) \
+    || defined(__BIG_ENDIAN__) || defined(BIG_ENDIAN)
+#      define DSFMT_BIG_ENDIAN 1
+#  endif
+#endif
+
+#if defined(DSFMT_BIG_ENDIAN) && defined(__amd64)
+#  undef DSFMT_BIG_ENDIAN
+#endif
+
+#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
+#  include <inttypes.h>
+#elif defined(_MSC_VER) || defined(__BORLANDC__)
+#  if !defined(DSFMT_UINT32_DEFINED) && !defined(SFMT_UINT32_DEFINED)
+typedef unsigned int uint32_t;
+typedef unsigned __int64 uint64_t;
+#    define UINT64_C(v) (v ## ui64)
+#    define DSFMT_UINT32_DEFINED
+#    if !defined(inline)
+#      define inline __inline
+#    endif
+#  endif
+#else
+#  include <inttypes.h>
+#  if !defined(inline)
+#    if defined(__GNUC__)
+#      define inline __inline__
+#    else
+#      define inline
+#    endif
+#  endif
+#endif
+
+#ifndef PRIu64
+#  if defined(_MSC_VER) || defined(__BORLANDC__)
+#    define PRIu64 "I64u"
+#    define PRIx64 "I64x"
+#  else
+#    define PRIu64 "llu"
+#    define PRIx64 "llx"
+#  endif
+#endif
+
+#ifndef UINT64_C
+#  define UINT64_C(v) (v ## ULL) 
+#endif
+
+/*------------------------------------------
+  128-bit SIMD like data type for standard C
+  ------------------------------------------*/
+#if defined(HAVE_ALTIVEC)
+#  if !defined(__APPLE__)
+#    include <altivec.h>
+#  endif
+/** 128-bit data structure */
+union W128_T {
+    vector unsigned int s;
+    uint64_t u[2];
+    uint32_t u32[4];
+    double d[2];
+};
+
+#elif defined(HAVE_SSE2)
+#  include <emmintrin.h>
+
+/** 128-bit data structure */
+union W128_T {
+    __m128i si;
+    __m128d sd;
+    uint64_t u[2];
+    uint32_t u32[4];
+    double d[2];
+};
+#else  /* standard C */
+/** 128-bit data structure */
+union W128_T {
+    uint64_t u[2];
+    uint32_t u32[4];
+    double d[2];
+};
+#endif
+
+/** 128-bit data type */
+typedef union W128_T w128_t;
+
+/** the 128-bit internal state array */
+struct DSFMT_T {
+    w128_t status[DSFMT_N + 1];
+    int idx;
+};
+typedef struct DSFMT_T dsfmt_t;
+
+/** dsfmt internal state vector */
+extern dsfmt_t dsfmt_global_data;
+/** dsfmt mexp for check */
+extern const int dsfmt_global_mexp;
+
+/* --- Ajout par Reza / Compilation C++ --- */ 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void dsfmt_gen_rand_all(dsfmt_t *dsfmt);
+void dsfmt_fill_array_open_close(dsfmt_t *dsfmt, double array[], int size);
+void dsfmt_fill_array_close_open(dsfmt_t *dsfmt, double array[], int size);
+void dsfmt_fill_array_open_open(dsfmt_t *dsfmt, double array[], int size);
+void dsfmt_fill_array_close1_open2(dsfmt_t *dsfmt, double array[], int size);
+void dsfmt_chk_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed, int mexp);
+void dsfmt_chk_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[],
+			     int key_length, int mexp);
+const char *dsfmt_get_idstring(void);
+int dsfmt_get_min_array_size(void);
+/* --- Ajout par Reza / Compilation C++ --- */ 
+#ifdef __cplusplus
+}
+#endif
+
+#if defined(__GNUC__)
+#  define DSFMT_PRE_INLINE inline static
+#  define DSFMT_PST_INLINE __attribute__((always_inline))
+#elif defined(_MSC_VER) && _MSC_VER >= 1200
+#  define DSFMT_PRE_INLINE __forceinline static
+#  define DSFMT_PST_INLINE
+#else
+#  define DSFMT_PRE_INLINE inline static
+#  define DSFMT_PST_INLINE
+#endif
+DSFMT_PRE_INLINE uint32_t dsfmt_genrand_uint32(dsfmt_t *dsfmt) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double dsfmt_genrand_close1_open2(dsfmt_t *dsfmt)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double dsfmt_genrand_close_open(dsfmt_t *dsfmt)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double dsfmt_genrand_open_close(dsfmt_t *dsfmt)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double dsfmt_genrand_open_open(dsfmt_t *dsfmt)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE uint32_t dsfmt_gv_genrand_uint32(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double dsfmt_gv_genrand_close1_open2(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double dsfmt_gv_genrand_close_open(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double dsfmt_gv_genrand_open_close(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double dsfmt_gv_genrand_open_open(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void dsfmt_gv_fill_array_open_close(double array[], int size)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void dsfmt_gv_fill_array_close_open(double array[], int size)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void dsfmt_gv_fill_array_open_open(double array[], int size)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void dsfmt_gv_fill_array_close1_open2(double array[], int size)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void dsfmt_gv_init_gen_rand(uint32_t seed) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void dsfmt_gv_init_by_array(uint32_t init_key[],
+					     int key_length) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void dsfmt_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void dsfmt_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[],
+					  int key_length) DSFMT_PST_INLINE;
+
+/**
+ * This function generates and returns unsigned 32-bit integer.
+ * This is slower than SFMT, only for convenience usage.
+ * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called
+ * before this function.
+ * @param dsfmt dsfmt internal state date
+ * @return double precision floating point pseudorandom number
+ */
+inline static uint32_t dsfmt_genrand_uint32(dsfmt_t *dsfmt) {
+    uint32_t r;
+    uint64_t *psfmt64 = &dsfmt->status[0].u[0];
+
+    if (dsfmt->idx >= DSFMT_N64) {
+	dsfmt_gen_rand_all(dsfmt);
+	dsfmt->idx = 0;
+    }
+    r = psfmt64[dsfmt->idx++] & 0xffffffffU;
+    return r;
+}
+
+/**
+ * This function generates and returns double precision pseudorandom
+ * number which distributes uniformly in the range [1, 2).  This is
+ * the primitive and faster than generating numbers in other ranges.
+ * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called
+ * before this function.
+ * @param dsfmt dsfmt internal state date
+ * @return double precision floating point pseudorandom number
+ */
+inline static double dsfmt_genrand_close1_open2(dsfmt_t *dsfmt) {
+    double r;
+    double *psfmt64 = &dsfmt->status[0].d[0];
+
+    if (dsfmt->idx >= DSFMT_N64) {
+	dsfmt_gen_rand_all(dsfmt);
+	dsfmt->idx = 0;
+    }
+    r = psfmt64[dsfmt->idx++];
+    return r;
+}
+
+/**
+ * This function generates and returns unsigned 32-bit integer.
+ * This is slower than SFMT, only for convenience usage.
+ * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called
+ * before this function.  This function uses \b global variables.
+ * @return double precision floating point pseudorandom number
+ */
+inline static uint32_t dsfmt_gv_genrand_uint32(void) {
+    return dsfmt_genrand_uint32(&dsfmt_global_data);
+}
+
+/**
+ * This function generates and returns double precision pseudorandom
+ * number which distributes uniformly in the range [1, 2).
+ * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called
+ * before this function. This function uses \b global variables.
+ * @return double precision floating point pseudorandom number
+ */
+inline static double dsfmt_gv_genrand_close1_open2(void) {
+    return dsfmt_genrand_close1_open2(&dsfmt_global_data);
+}
+
+/**
+ * This function generates and returns double precision pseudorandom
+ * number which distributes uniformly in the range [0, 1).
+ * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called
+ * before this function.
+ * @param dsfmt dsfmt internal state date
+ * @return double precision floating point pseudorandom number
+ */
+inline static double dsfmt_genrand_close_open(dsfmt_t *dsfmt) {
+    return dsfmt_genrand_close1_open2(dsfmt) - 1.0;
+}
+
+/**
+ * This function generates and returns double precision pseudorandom
+ * number which distributes uniformly in the range [0, 1).
+ * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called
+ * before this function. This function uses \b global variables.
+ * @return double precision floating point pseudorandom number
+ */
+inline static double dsfmt_gv_genrand_close_open(void) {
+    return dsfmt_gv_genrand_close1_open2() - 1.0;
+}
+
+/**
+ * This function generates and returns double precision pseudorandom
+ * number which distributes uniformly in the range (0, 1].
+ * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called
+ * before this function. 
+ * @param dsfmt dsfmt internal state date
+ * @return double precision floating point pseudorandom number
+ */
+inline static double dsfmt_genrand_open_close(dsfmt_t *dsfmt) {
+    return 2.0 - dsfmt_genrand_close1_open2(dsfmt);
+}
+
+/**
+ * This function generates and returns double precision pseudorandom
+ * number which distributes uniformly in the range (0, 1].
+ * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called
+ * before this function. This function uses \b global variables.
+ * @return double precision floating point pseudorandom number
+ */
+inline static double dsfmt_gv_genrand_open_close(void) {
+    return 2.0 - dsfmt_gv_genrand_close1_open2();
+}
+
+/**
+ * This function generates and returns double precision pseudorandom
+ * number which distributes uniformly in the range (0, 1).
+ * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called
+ * before this function.
+ * @param dsfmt dsfmt internal state date
+ * @return double precision floating point pseudorandom number
+ */
+inline static double dsfmt_genrand_open_open(dsfmt_t *dsfmt) {
+    double *dsfmt64 = &dsfmt->status[0].d[0];
+    union {
+	double d;
+	uint64_t u;
+    } r;
+
+    if (dsfmt->idx >= DSFMT_N64) {
+	dsfmt_gen_rand_all(dsfmt);
+	dsfmt->idx = 0;
+    }
+    r.d = dsfmt64[dsfmt->idx++];
+    r.u |= 1;
+    return r.d - 1.0;
+}
+
+/**
+ * This function generates and returns double precision pseudorandom
+ * number which distributes uniformly in the range (0, 1).
+ * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called
+ * before this function. This function uses \b global variables.
+ * @return double precision floating point pseudorandom number
+ */
+inline static double dsfmt_gv_genrand_open_open(void) {
+    return dsfmt_genrand_open_open(&dsfmt_global_data);
+}
+
+/**
+ * This function generates double precision floating point
+ * pseudorandom numbers which distribute in the range [1, 2) to the
+ * specified array[] by one call. This function is the same as
+ * dsfmt_fill_array_close1_open2() except that this function uses
+ * \b global variables.
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa dsfmt_fill_array_close1_open2()
+ */
+inline static void dsfmt_gv_fill_array_close1_open2(double array[], int size) {
+    dsfmt_fill_array_close1_open2(&dsfmt_global_data, array, size);
+}
+
+/**
+ * This function generates double precision floating point
+ * pseudorandom numbers which distribute in the range (0, 1] to the
+ * specified array[] by one call. This function is the same as
+ * dsfmt_gv_fill_array_close1_open2() except the distribution range.
+ * This function uses \b global variables.
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa dsfmt_fill_array_close1_open2() and \sa
+ * dsfmt_gv_fill_array_close1_open2()
+ */
+inline static void dsfmt_gv_fill_array_open_close(double array[], int size) {
+    dsfmt_fill_array_open_close(&dsfmt_global_data, array, size);
+}
+
+/**
+ * This function generates double precision floating point
+ * pseudorandom numbers which distribute in the range [0, 1) to the
+ * specified array[] by one call. This function is the same as
+ * dsfmt_gv_fill_array_close1_open2() except the distribution range.
+ * This function uses \b global variables.
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa dsfmt_fill_array_close1_open2() \sa
+ * dsfmt_gv_fill_array_close1_open2()
+ */
+inline static void dsfmt_gv_fill_array_close_open(double array[], int size) {
+    dsfmt_fill_array_close_open(&dsfmt_global_data, array, size);
+}
+
+/**
+ * This function generates double precision floating point
+ * pseudorandom numbers which distribute in the range (0, 1) to the
+ * specified array[] by one call. This function is the same as
+ * dsfmt_gv_fill_array_close1_open2() except the distribution range.
+ * This function uses \b global variables.
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa dsfmt_fill_array_close1_open2() \sa
+ * dsfmt_gv_fill_array_close1_open2()
+ */
+inline static void dsfmt_gv_fill_array_open_open(double array[], int size) {
+    dsfmt_fill_array_open_open(&dsfmt_global_data, array, size);
+}
+
+/**
+ * This function initializes the internal state array with a 32-bit
+ * integer seed.
+ * @param dsfmt dsfmt state vector.
+ * @param seed a 32-bit integer used as the seed.
+ */
+inline static void dsfmt_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed) {
+    dsfmt_chk_init_gen_rand(dsfmt, seed, DSFMT_MEXP);
+}
+
+/**
+ * This function initializes the internal state array with a 32-bit
+ * integer seed. This function uses \b global variables.
+ * @param seed a 32-bit integer used as the seed.
+ * see also \sa dsfmt_init_gen_rand()
+ */
+inline static void dsfmt_gv_init_gen_rand(uint32_t seed) {
+    dsfmt_init_gen_rand(&dsfmt_global_data, seed);
+}
+
+/**
+ * This function initializes the internal state array,
+ * with an array of 32-bit integers used as the seeds.
+ * @param dsfmt dsfmt state vector
+ * @param init_key the array of 32-bit integers, used as a seed.
+ * @param key_length the length of init_key.
+ */
+inline static void dsfmt_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[],
+				       int key_length) {
+    dsfmt_chk_init_by_array(dsfmt, init_key, key_length, DSFMT_MEXP);
+}
+
+/**
+ * This function initializes the internal state array,
+ * with an array of 32-bit integers used as the seeds.
+ * This function uses \b global variables.
+ * @param init_key the array of 32-bit integers, used as a seed.
+ * @param key_length the length of init_key.
+ * see also \sa dsfmt_init_by_array()
+ */
+inline static void dsfmt_gv_init_by_array(uint32_t init_key[], int key_length) {
+    dsfmt_init_by_array(&dsfmt_global_data, init_key, key_length);
+}
+
+#if !defined(DSFMT_DO_NOT_USE_OLD_NAMES)
+DSFMT_PRE_INLINE const char *get_idstring(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE int get_min_array_size(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void init_gen_rand(uint32_t seed) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void init_by_array(uint32_t init_key[], int key_length)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double genrand_close1_open2(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double genrand_close_open(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double genrand_open_close(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE double genrand_open_open(void) DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void fill_array_open_close(double array[], int size)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void fill_array_close_open(double array[], int size)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void fill_array_open_open(double array[], int size)
+    DSFMT_PST_INLINE;
+DSFMT_PRE_INLINE void fill_array_close1_open2(double array[], int size)
+    DSFMT_PST_INLINE;
+
+/**
+ * This function is just the same as dsfmt_get_idstring().
+ * @return id string.
+ * see also \sa dsfmt_get_idstring()
+ */
+inline static const char *get_idstring(void) {
+    return dsfmt_get_idstring();
+}
+
+/**
+ * This function is just the same as dsfmt_get_min_array_size().
+ * @return minimum size of array used for fill_array functions.
+ * see also \sa dsfmt_get_min_array_size()
+ */
+inline static int get_min_array_size(void) {
+    return dsfmt_get_min_array_size();
+}
+
+/**
+ * This function is just the same as dsfmt_gv_init_gen_rand().
+ * @param seed a 32-bit integer used as the seed.
+ * see also \sa dsfmt_gv_init_gen_rand(), \sa dsfmt_init_gen_rand().
+ */
+inline static void init_gen_rand(uint32_t seed) {
+    dsfmt_gv_init_gen_rand(seed);
+}
+
+/**
+ * This function is just the same as dsfmt_gv_init_by_array().
+ * @param init_key the array of 32-bit integers, used as a seed.
+ * @param key_length the length of init_key.
+ * see also \sa dsfmt_gv_init_by_array(), \sa dsfmt_init_by_array().
+ */
+inline static void init_by_array(uint32_t init_key[], int key_length) {
+    dsfmt_gv_init_by_array(init_key, key_length);
+}
+
+/**
+ * This function is just the same as dsfmt_gv_genrand_close1_open2().
+ * @return double precision floating point number.
+ * see also \sa dsfmt_genrand_close1_open2() \sa
+ * dsfmt_gv_genrand_close1_open2()
+ */
+inline static double genrand_close1_open2(void) {
+    return dsfmt_gv_genrand_close1_open2();
+}
+
+/**
+ * This function is just the same as dsfmt_gv_genrand_close_open().
+ * @return double precision floating point number.
+ * see also \sa dsfmt_genrand_close_open() \sa
+ * dsfmt_gv_genrand_close_open()
+ */
+inline static double genrand_close_open(void) {
+    return dsfmt_gv_genrand_close_open();
+}
+
+/**
+ * This function is just the same as dsfmt_gv_genrand_open_close().
+ * @return double precision floating point number.
+ * see also \sa dsfmt_genrand_open_close() \sa
+ * dsfmt_gv_genrand_open_close()
+ */
+inline static double genrand_open_close(void) {
+    return dsfmt_gv_genrand_open_close();
+}
+
+/**
+ * This function is just the same as dsfmt_gv_genrand_open_open().
+ * @return double precision floating point number.
+ * see also \sa dsfmt_genrand_open_open() \sa
+ * dsfmt_gv_genrand_open_open()
+ */
+inline static double genrand_open_open(void) {
+    return dsfmt_gv_genrand_open_open();
+}
+
+/**
+ * This function is juset the same as dsfmt_gv_fill_array_open_close().
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa dsfmt_gv_fill_array_open_close(), \sa
+ * dsfmt_fill_array_close1_open2(), \sa
+ * dsfmt_gv_fill_array_close1_open2()
+ */
+inline static void fill_array_open_close(double array[], int size) {
+    dsfmt_gv_fill_array_open_close(array, size);
+}
+
+/**
+ * This function is juset the same as dsfmt_gv_fill_array_close_open().
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa dsfmt_gv_fill_array_close_open(), \sa
+ * dsfmt_fill_array_close1_open2(), \sa
+ * dsfmt_gv_fill_array_close1_open2()
+ */
+inline static void fill_array_close_open(double array[], int size) {
+    dsfmt_gv_fill_array_close_open(array, size);
+}
+
+/**
+ * This function is juset the same as dsfmt_gv_fill_array_open_open().
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa dsfmt_gv_fill_array_open_open(), \sa
+ * dsfmt_fill_array_close1_open2(), \sa
+ * dsfmt_gv_fill_array_close1_open2()
+ */
+inline static void fill_array_open_open(double array[], int size) {
+    dsfmt_gv_fill_array_open_open(array, size);
+}
+
+/**
+ * This function is juset the same as dsfmt_gv_fill_array_close1_open2().
+ * @param array an array where pseudorandom numbers are filled
+ * by this function.
+ * @param size the number of pseudorandom numbers to be generated.
+ * see also \sa dsfmt_fill_array_close1_open2(), \sa
+ * dsfmt_gv_fill_array_close1_open2()
+ */
+inline static void fill_array_close1_open2(double array[], int size) {
+    dsfmt_gv_fill_array_close1_open2(array, size);
+}
+#endif /* DSFMT_DO_NOT_USE_OLD_NAMES */
+
+#endif /* DSFMT_H */
Index: trunk/SophyaLib/BaseTools/dsfmtflags.h
===================================================================
--- trunk/SophyaLib/BaseTools/dsfmtflags.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/dsfmtflags.h	(revision 3602)
@@ -0,0 +1,19 @@
+#ifndef  DSFMTWRAP_H_SEEN 
+#define  DSFMTWRAP_H_SEEN 
+
+// Fichier include d'enrobage pour flags et definition du generateur 
+// SIMD-oriented Fast Mersenne Twister (SFMT).
+// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+
+// Definir les flags en fonction des architectures ...
+// Choix de la periode du generateur 
+//#define  DSFMT_MEXP  521
+#define  DSFMT_MEXP  19937
+
+// Si on veut generer les instructions SSE2, : 
+// Streaming SIMD Extensions 2, is one of the IA-32 SIMD (Single Instruction, Multiple Data) 
+#if defined(MACH_INTEL) || defined(MACH_AMD) 
+#define  HAVE_SSE2 1
+#endif
+
+#endif
Index: trunk/SophyaLib/BaseTools/objlist.list
===================================================================
--- trunk/SophyaLib/BaseTools/objlist.list	(revision 3601)
+++ trunk/SophyaLib/BaseTools/objlist.list	(revision 3602)
@@ -1,3 +1,4 @@
 datatype.o
+dSFMT.o
 dvlist.o
 fiondblock.o
@@ -9,4 +10,7 @@
 ppfbinstream.o
 ppftpointerio.o
+randfmt.o
+randinterf.o
+randr48.o
 rawstream.o
 sophyainit.o
Index: trunk/SophyaLib/BaseTools/randfmt.cc
===================================================================
--- trunk/SophyaLib/BaseTools/randfmt.cc	(revision 3602)
+++ trunk/SophyaLib/BaseTools/randfmt.cc	(revision 3602)
@@ -0,0 +1,100 @@
+#include <math.h> 
+#include <stdlib.h>
+
+#include "randfmt.h"
+
+namespace SOPHYA {
+
+//----------------------------------------------------------------
+// Implementation d'un generateur aleatoire utilisant le code 
+// SIMD-oriented Fast Mersenne Twister (SFMT).
+// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+
+
+FMTRandGen::FMTRandGen(uint_4 seed)
+{
+  dsfmt_init_gen_rand(&dsfmt_,seed);
+}
+
+FMTRandGen::~FMTRandGen()
+{
+}
+
+void FMTRandGen::SetSeed(uint_4 seed)
+{
+  dsfmt_init_gen_rand(&dsfmt_,seed);
+}
+
+void FMTRandGen::SetSeed(vector<uint_4> seed)
+{
+  if(seed.size()<=0) return;
+  int key_length = (int)seed.size();
+  uint_4* init_key = new uint_4[key_length];
+  for(int i=0;i<key_length;i++) init_key[i] = seed[i];
+
+  dsfmt_init_by_array(&dsfmt_,init_key,key_length);
+
+  delete [] init_key;
+}
+
+r_8 FMTRandGen::Next()
+{
+  return dsfmt_genrand_close_open(&dsfmt_);
+}
+
+void FMTRandGen::AutoInit(int lp)
+{
+  vector<uint_2> seed;
+  GenerateSeedVector(1,seed,lp);
+  vector<uint_4> s;
+  uint_4 s4;
+  s4 = seed[0] | (seed[1]<<16);
+  s.push_back(s4);
+  s4 = seed[2] | (seed[3]<<16);
+  s.push_back(s4);
+  SetSeed(s);
+}
+
+//----------------------------------------------------------
+// Classe pour la gestion de persistance
+// ObjFileIO<FMTRandGen>
+//----------------------------------------------------------
+
+/* --Methode-- */
+DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
+void ObjFileIO<FMTRandGen>::WriteSelf(POutPersist& s) const
+{
+  if (dobj == NULL)
+    throw NullPtrError("ObjFileIO<FMTRandGen>::WriteSelf() dobj=NULL");
+  s.Put((uint_4)DSFMT_MEXP);
+  s.Put(dobj->dsfmt_.idx);
+  for(int i=0;i<DSFMT_N+1;i++) 
+    for(int j=0;j<2;j++) s.Put(dobj->dsfmt_.status[i].u[j]);
+  return;
+}
+
+/* --Methode-- */
+DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
+void ObjFileIO<FMTRandGen>::ReadSelf(PInPersist& s)
+{
+  uint_4 mexp;
+  s.Get(mexp);
+  if(mexp != DSFMT_MEXP)
+    throw SzMismatchError("ObjFileIO<FMTRandGen>::WriteSelf() wrong DSFMT_MEXP");
+  if(dobj == NULL) dobj = new FMTRandGen();
+  s.Get(dobj->dsfmt_.idx);
+  for(int i=0;i<DSFMT_N+1;i++) 
+    for(int j=0;j<2;j++) s.Get(dobj->dsfmt_.status[i].u[j]);
+  return;
+}
+
+// ---------------------------------------------------------
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template ObjFileIO<FMTRandGen>
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class ObjFileIO<FMTRandGen>;
+#endif
+
+}  /* namespace SOPHYA */
Index: trunk/SophyaLib/BaseTools/randfmt.h
===================================================================
--- trunk/SophyaLib/BaseTools/randfmt.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/randfmt.h	(revision 3602)
@@ -0,0 +1,52 @@
+#ifndef FMTRANDGEN_H_SEEN
+#define FMTRANDGEN_H_SEEN
+
+// Classes FMTRandGen
+// Generateur aleatoire utilisant le code dSFMT 
+//
+// R. Ansari          (C) UPS+LAL IN2P3/CNRS  2009
+// C. Magneville      (C) DAPNIA/SPP  CEA     2009
+
+
+#include "machdefs.h"
+#include "objfio.h"
+
+#include "dsfmtflags.h"
+#include "dSFMT.h"
+#include "randinterf.h"
+
+#include <vector>
+
+namespace SOPHYA {
+
+//! Pseudorandom number generator class using dSFMT code 
+class FMTRandGen : public RandomGeneratorInterface {
+
+ public:
+  FMTRandGen(uint_4 seed=12345);
+  virtual ~FMTRandGen();
+
+  virtual void SetSeed(uint_4 seed);
+  virtual void SetSeed(vector<uint_4> seed);
+  virtual void AutoInit(int lp=0);
+
+  friend class ObjFileIO<FMTRandGen> ;
+
+ protected:
+  virtual r_8 Next();
+  dsfmt_t dsfmt_;
+
+};
+
+// Classe pour la gestion de persistance PPF :  ObjFileIO<FMTRandGen>
+
+/*! Writes the random generator object state in the POutPersist stream \b os */
+inline POutPersist& operator << (POutPersist& os, FMTRandGen & obj)
+{ ObjFileIO<FMTRandGen> fio(&obj);  fio.Write(os);  return(os); }
+/*! Reads the random generator object state from the PInPersist stream \b is */
+inline PInPersist& operator >> (PInPersist& is, FMTRandGen & obj)
+{ ObjFileIO<FMTRandGen> fio(&obj); is.SkipToNextObject(); fio.Read(is); return(is); }
+
+} /* namespace SOPHYA */
+
+#endif 
Index: trunk/SophyaLib/BaseTools/randinterf.cc
===================================================================
--- trunk/SophyaLib/BaseTools/randinterf.cc	(revision 3602)
+++ trunk/SophyaLib/BaseTools/randinterf.cc	(revision 3602)
@@ -0,0 +1,752 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <math.h> 
+#include <stdlib.h> 
+#include <sys/time.h>
+#include <time.h>
+#include <iostream>
+#include "pexceptions.h"
+
+#include "randinterf.h"
+
+namespace SOPHYA {
+
+//-------------------------------------------------------------------------------
+// ------ Definition d'interface des classes de generateurs de nombres aleatoires
+
+RandomGeneratorInterface::RandomGeneratorInterface()
+  : usegaussian_(0), usepoisson_(0), useexpo_(0)
+{
+}
+ 
+
+RandomGeneratorInterface::~RandomGeneratorInterface(void) 
+{
+  // rien a faire 
+}
+
+void RandomGeneratorInterface::UseGaussian(uint_2 type)
+{
+  usegaussian_ = type;
+}
+
+void RandomGeneratorInterface::UseExpo(uint_2 type)
+{
+  useexpo_ = type;
+}
+
+void RandomGeneratorInterface::UsePoisson(uint_2 type)
+{
+  usepoisson_ = type;
+}
+
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+
+r_8 RandomGeneratorInterface::Next()
+{
+  printf("RandomGeneratorInterface::Next(): undefined code !!!\n");
+  throw MathExc("RandomGeneratorInterface::Next(): undefined code !!!");
+}
+
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+void RandomGeneratorInterface::GenerateSeedVector(int nseed,vector<uint_2>& seed,int lp)
+// renvoie un vecteur de nseed+2 entiers 32 bits
+// [0 - 2] = codage sur 48 bits du nombre (melange) de microsec depuis l'origine
+// [3 -> 3+ngene-1] = entiers aleatoires (poor man generator)
+//
+// L'initialiseur est donne par un codage du nombre de millisecondes
+// ecoulees depuis le 0 heure le 1er Janvier 1970 UTC (cf gettimeofday).
+// Seuls les 48 bits de poids faible sont retenus.
+// Un melange des bits est ensuite effectue pour que les 3 nombres
+// (unsigned short) d'initialisation ne soient pas trop semblables.
+// Le nombre le plus grand que l'on peut mettre
+// dans un entier unsigned de N bits est: 2^N-1
+// 48 bits -> 2^48-1 = 281474976710655 musec = 3257.8j = 8.9y
+//         -> meme initialisation tous les 8.9 ans a 1 microsec pres !
+{
+  if(lp>0) cout<<"RandomGeneratorInterface::GenerateSeedVector: nseed="<<nseed<<endl;
+
+  // ---
+  // --- les deux premiers mots remplis avec le temps
+  // ---
+  // On recupere le temps ecoule depuis l'origine code en sec+musec
+  struct timeval now;
+  gettimeofday(&now,0);
+  // Calcul du temps ecoule depuis l'origine en microsecondes
+  uint_8 tmicro70 = (uint_8)now.tv_sec*(uint_8)1000000 + (uint_8)now.tv_usec;
+  if(lp>1) cout<<"."<<now.tv_sec<<" sec + "<<now.tv_usec<<" musec = "<<tmicro70<<" musec"<<endl;
+  // Remplissage du tableau de 48 bits
+  uint_2 b[48]; uint_8 tdum = tmicro70;
+  for(int ip=0;ip<48;ip++) {b[ip] = tdum&1; tdum = (tdum>>1);}
+  if(lp>2) {
+    cout<<"..b= ";
+    for(int ip=47;ip>=0;ip--) {cout<<b[ip]; if(ip%32==0 || ip%16==0) cout<<" ";}
+    cout<<endl;
+  }
+  // Melange des bits qui varient vite (poids faible, microsec)
+  //   avec ceux variant lentement (poids fort, sec)
+  for(int ip=0;ip<16;ip++) {
+    if(ip%3==1) swap(b[ip],b[32+ip]);
+    else if(ip%3==2) swap(b[ip],b[16-ip]);
+  }
+  if(lp>2) {
+    cout<<"..b= ";
+    for(int ip=47;ip>=0;ip--) {cout<<b[ip]; if(ip%32==0 || ip%16==0) cout<<" ";}
+    cout<<endl;
+  }
+  // Remplissage
+  seed.resize(0);
+  for(int i=0;i<3;i++) {
+    seed.push_back(0);
+    uint_2 w = 1;
+    for(int ip=0;ip<16;ip++) {seed[i] += w*b[i*16+ip]; w *= 2;}
+  }
+  if(lp>0) cout<<"seed(time): "<<seed[0]<<" "<<seed[1]<<" "<<seed[2]<<endl;
+
+  // ---
+  // --- generation des nombres aleatoires complementaires (poor man generator)
+  // ---
+  //----------------------------------------------------------------------------//
+  // Ran088: L'Ecuyer's 1996 three-component Tausworthe generator "taus88"
+  // Returns an integer random number uniformly distributed within [0,4294967295]
+  // The period length is approximately 2^88 (which is 3*10^26). 
+  // This generator is very fast and passes all standard statistical tests.
+  // Reference:
+  //   (1) P. L'Ecuyer, Maximally equidistributed combined Tausworthe generators,
+  //       Mathematics of Computation, 65, 203-213 (1996), see Figure 4.
+  //   (2) recommended in:
+  //       P. L'Ecuyer, Random number generation, chapter 4 of the
+  //       Handbook on Simulation, Ed. Jerry Banks, Wiley, 1997.
+  //----------------------------------------------------------------------------//
+  if(nseed<=0) return;
+  // initialize seeds using the given seed value taking care of
+  // the requirements. The constants below are arbitrary otherwise
+  uint_4 seed0 = uint_4(tmicro70&0xFFFFFFFFULL);
+  if(lp>2) cout<<"seed0(time): "<<seed0<<endl;
+  uint_4 state_s1, state_s2, state_s3;
+  state_s1 = 1243598713U ^ seed0; if (state_s1 <  2) state_s1 = 1243598713U;
+  state_s2 = 3093459404U ^ seed0; if (state_s2 <  8) state_s2 = 3093459404U;
+  state_s3 = 1821928721U ^ seed0; if (state_s3 < 16) state_s3 = 1821928721U;
+  int nfill = 0, ico=0;	
+  while(nfill<nseed) {
+    uint_4 s1 = state_s1, s2 = state_s2, s3 = state_s3;
+    // generate a random 32 bit number
+    s1 = ((s1 &  -2) << 12) ^ (((s1 << 13) ^  s1) >> 19);
+    s2 = ((s2 &  -8) <<  4) ^ (((s2 <<  2) ^  s2) >> 25);
+    s3 = ((s3 & -16) << 17) ^ (((s3 <<  3) ^  s3) >> 11);
+    state_s1 = s1; state_s2 = s2; state_s3 = s3;
+    // le nombre aleatoire sur 32 bits est: s1^s2^s3
+    if(ico<15) {ico++; continue;}  // des tirages blancs
+    uint_2 s = uint_2( (s1^s2^s3)&0xFFFFU );
+    seed.push_back(s);
+    if(lp>0) cout<<"seed(t88): "<<seed[3+nfill]<<endl;
+    nfill++;
+  }
+
+}
+
+
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+
+//--- Generation de nombre aleatoires suivant une distribution gaussienne 
+r_8 RandomGeneratorInterface::GaussianBoxMuller() 
+{
+  r_8 A=Next(); 
+  while (A==0.) A=Next();
+  return sqrt(-2.*log(A))*cos(2.*M_PI*Next());
+}
+
+//-------------------------------------------
+// Adapte de ranlib float snorm() 
+// http://orion.math.iastate.edu/burkardt/c_src/ranlib/ranlib.c
+/*
+**********************************************************************
+     (STANDARD-)  N O R M A L  DISTRIBUTION
+**********************************************************************
+                                                                      
+     FOR DETAILS SEE:                                                 
+                                                                      
+               AHRENS, J.H. AND DIETER, U.                            
+               EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             
+               SAMPLING FROM THE NORMAL DISTRIBUTION.                 
+               MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          
+                                                                      
+     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  
+     (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  
+                                                                      
+     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
+     SUNIF.  The argument IR thus goes away.                          
+                                                                      
+**********************************************************************
+     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
+     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
+*/
+static double a_snorm[32] = {
+    0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
+    0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
+    0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
+    1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
+    1.862732,2.153875
+};
+static double d_snorm[31] = {
+    0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
+    0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
+    0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
+    0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
+};
+static float t_snorm[31] = {
+    7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
+    1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
+    2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
+    4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
+    9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
+};
+static float h_snorm[31] = {
+    3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
+    4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
+    4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
+    5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
+    8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
+};
+r_8 RandomGeneratorInterface::GaussianSNorm()
+{
+long i;
+double snorm,u,s,ustar,aa,w,y,tt;
+    u = Next();
+    s = 0.0;
+    if(u > 0.5) s = 1.0;
+    u += (u-s);
+    u = 32.0*u;
+    i = (long) (u);
+    if(i == 32) i = 31;
+    if(i == 0) goto S100;
+/*
+                                START CENTER
+*/
+    ustar = u-(double)i;
+    aa = *(a_snorm+i-1);
+S40:
+    if(ustar <= *(t_snorm+i-1)) goto S60;
+    w = (ustar-*(t_snorm+i-1))**(h_snorm+i-1);
+S50:
+/*
+                                EXIT   (BOTH CASES)
+*/
+    y = aa+w;
+    snorm = y;
+    if(s == 1.0) snorm = -y;
+    return snorm;
+S60:
+/*
+                                CENTER CONTINUED
+*/
+    u = Next();
+    w = u*(*(a_snorm+i)-aa);
+    tt = (0.5*w+aa)*w;
+    goto S80;
+S70:
+    tt = u;
+    ustar = Next();
+S80:
+    if(ustar > tt) goto S50;
+    u = Next();
+    if(ustar >= u) goto S70;
+    ustar = Next();
+    goto S40;
+S100:
+/*
+                                START TAIL
+*/
+    i = 6;
+    aa = *(a_snorm+31);
+    goto S120;
+S110:
+    aa += *(d_snorm+i-1);
+    i += 1;
+S120:
+    u += u;
+    if(u < 1.0) goto S110;
+    u -= 1.0;
+S140:
+    w = u**(d_snorm+i-1);
+    tt = (0.5*w+aa)*w;
+    goto S160;
+S150:
+    tt = u;
+S160:
+    ustar = Next();
+    if(ustar > tt) goto S50;
+    u = Next();
+    if(ustar >= u) goto S150;
+    u = Next();
+    goto S140;
+}
+
+r_8 RandomGeneratorInterface::GaussianPolarBoxMuller()
+{
+double x1,x2,w;
+do {
+   x1 = 2.0 * Next() - 1.0;
+   x2 = 2.0 * Next() - 1.0;
+   w = x1 * x1 + x2 * x2;
+   } while ( w >= 1.0 || w==0. );
+return x1 * sqrt(-2.0*log(w)/w);
+}
+
+static double s2se_RatioUnif=sqrt(2./M_E) , epm135_RatioUnif=exp(-1.35) , ep1q_RatioUnif=exp(1./4.);
+r_8 RandomGeneratorInterface::GaussianRatioUnif()
+{
+double u,v,x;
+while(true) {
+  do {u = Next();} while ( u == 0. );
+  v = (2.0*Next()-1.0)*s2se_RatioUnif;
+  x = v/u;
+  if(x*x <= 5.0-4.0*ep1q_RatioUnif*u) break;
+  if(x*x<4.0*epm135_RatioUnif/u+1.4)
+    if(v*v<-4.0*u*u*log(u)) break;
+}
+return x;
+}
+
+r_8 RandomGeneratorInterface::GaussianLevaRatioUnif()
+{
+double u,v,x,y,q;
+ do {
+   u = 1.-Next();  // in ]0,1]
+   v = Next()-0.5;  // in [-0.5, 0.5[
+   v *= 1.7156;
+   x = u - 0.449871;
+   y = ((v<0)?-v:v) + 0.386595;
+   q = x*x + y*(0.19600*y - 0.25472*x);
+ } while( q>=0.27597 && (q>0.27846  || v*v>-4.0*u*u*log(u)) );
+ return v/u;
+}
+
+r_8 RandomGeneratorInterface::GaussianTail(double s)
+{
+  /* Returns a gaussian random variable larger than a
+   * This implementation does one-sided upper-tailed deviates.
+   */
+
+  if (s < 1)
+    {
+      /* For small s, use a direct rejection method. The limit s < 1
+         can be adjusted to optimise the overall efficiency */
+      double x;
+      do
+        {
+          x = Gaussian();
+        }
+      while (x < s);
+      return x;
+    }
+  else
+    {
+      /* Use the "supertail" deviates from the last two steps
+       * of Marsaglia's rectangle-wedge-tail method, as described
+       * in Knuth, v2, 3rd ed, pp 123-128.  (See also exercise 11, p139,
+       * and the solution, p586.)
+       */
+      double u, v, x;
+      do
+        {
+          u = Next();
+          do
+            {
+              v = Next();
+            }
+          while (v == 0.0);
+          x = sqrt (s * s - 2 * log (v));
+        }
+      while (x * u > s);
+      return x;
+    }
+}
+
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+
+//--- Generation de nombre aleatoires suivant une distribution de Poisson 
+uint_8 RandomGeneratorInterface::PoissonSimple(double mu,double mumax)
+{
+  double pp,ppi,x;
+
+  if((mumax>0.)&&(mu>=mumax)) {
+    pp = sqrt(mu);
+    while( (x=pp*Gaussian()) < -mu );
+    return (uint_8)(mu+x+0.5);
+  }
+  else {
+    uint_8 n;
+    ppi = pp = exp(-mu);
+    x = Next();
+    n = 0;
+    while (x > ppi) {
+      n++;
+      pp = mu*pp/(double)n;
+      ppi += pp;
+    }
+    return n;
+  }
+  return 0;  // pas necessaire ?
+}
+
+
+static double a0_poiahr = -0.5;
+static double a1_poiahr = 0.3333333;
+static double a2_poiahr = -0.2500068;
+static double a3_poiahr = 0.2000118;
+static double a4_poiahr = -0.1661269;
+static double a5_poiahr = 0.1421878;
+static double a6_poiahr = -0.1384794;
+static double a7_poiahr = 0.125006;
+static double fact_poiahr[10] = {
+    1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0};
+uint_8 RandomGeneratorInterface::PoissonAhrens(double mu)
+/*
+**********************************************************************
+     long ignpoi(float mu)
+                    GENerate POIsson random deviate
+                              Function
+     Generates a single random deviate from a Poisson
+     distribution with mean AV.
+                              Arguments
+     av --> The mean of the Poisson distribution from which
+            a random deviate is to be generated.
+     genexp <-- The random deviate.
+                              Method
+     Renames KPOIS from TOMS as slightly modified by BWB to use RANF
+     instead of SUNIF.
+     For details see:
+               Ahrens, J.H. and Dieter, U.
+               Computer Generation of Poisson Deviates
+               From Modified Normal Distributions.
+               ACM Trans. Math. Software, 8, 2
+               (June 1982),163-179
+**********************************************************************
+**********************************************************************
+                                                                      
+                                                                      
+     P O I S S O N  DISTRIBUTION                                      
+                                                                      
+                                                                      
+**********************************************************************
+**********************************************************************
+                                                                      
+     FOR DETAILS SEE:                                                 
+                                                                      
+               AHRENS, J.H. AND DIETER, U.                            
+               COMPUTER GENERATION OF POISSON DEVIATES                
+               FROM MODIFIED NORMAL DISTRIBUTIONS.                    
+               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. 
+                                                                      
+     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)  
+                                                                      
+**********************************************************************
+      INTEGER FUNCTION IGNPOI(IR,MU)
+     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
+             MU=MEAN MU OF THE POISSON DISTRIBUTION
+     OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
+     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
+     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
+     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
+     SEPARATION OF CASES A AND B
+*/
+{
+uint_8 long ignpoi,j,k,kflag,l,m;
+double b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g,omega,p,p0,px,py,q,s,
+    t,u,v,x,xx,pp[35];
+
+    if(mu < 10.0) goto S120;
+/*
+     C A S E  A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
+*/
+    s = sqrt(mu);
+    d = 6.0*mu*mu;
+/*
+             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
+             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
+             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
+*/
+    l = (uint_8) (mu-1.1484);
+/*
+     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
+*/
+    g = mu+s*Gaussian();
+    if(g < 0.0) goto S20;
+    ignpoi = (uint_8) (g);
+/*
+     STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
+*/
+    if(ignpoi >= l) return ignpoi;
+/*
+     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
+*/
+    fk = (double)ignpoi;
+    difmuk = mu-fk;
+    u = Next();
+    if(d*u >= difmuk*difmuk*difmuk) return ignpoi;
+S20:
+/*
+     STEP P. PREPARATIONS FOR STEPS Q AND H.
+             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
+             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
+             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
+             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
+             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
+*/
+    omega = 0.3989423/s;
+    b1 = 4.166667E-2/mu;
+    b2 = 0.3*b1*b1;
+    c3 = 0.1428571*b1*b2;
+    c2 = b2-15.0*c3;
+    c1 = b1-6.0*b2+45.0*c3;
+    c0 = 1.0-b1+3.0*b2-15.0*c3;
+    c = 0.1069/mu;
+    if(g < 0.0) goto S50;
+/*
+             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
+*/
+    kflag = 0;
+    goto S70;
+S40:
+/*
+     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
+*/
+    if(fy-u*fy <= py*exp(px-fx)) return ignpoi;
+S50:
+/*
+     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
+             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
+             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
+*/
+    e = Expo();
+    u = Next();
+    u += (u-1.0);
+    //t = 1.8+fsign(e,u);
+    t = 1.8 + (((u>0. && e<0.) || (u<0. && e>0.))?-e:e);
+    if(t <= -0.6744) goto S50;
+    ignpoi = (uint_8) (mu+s*t);
+    fk = (double)ignpoi;
+    difmuk = mu-fk;
+/*
+             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
+*/
+    kflag = 1;
+    goto S70;
+S60:
+/*
+     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
+*/
+    if(c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50;
+    return ignpoi;
+S70:
+/*
+     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
+             CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
+*/
+    if(ignpoi >= 10) goto S80;
+    px = -mu;
+    py = pow(mu,(double)ignpoi)/ *(fact_poiahr+ignpoi);
+    goto S110;
+S80:
+/*
+             CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
+             A0-A7 FOR ACCURACY WHEN ADVISABLE
+             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
+*/
+    del = 8.333333E-2/fk;
+    del -= (4.8*del*del*del);
+    v = difmuk/fk;
+    if(fabs(v) <= 0.25) goto S90;
+    px = fk*log(1.0+v)-difmuk-del;
+    goto S100;
+S90:
+    px = fk*v*v*(((((((a7_poiahr*v+a6_poiahr)*v+a5_poiahr)*v+a4_poiahr)*v+a3_poiahr)*v+a2_poiahr)*v+a1_poiahr)*v+a0_poiahr)-del;
+S100:
+    py = 0.3989423/sqrt(fk);
+S110:
+    x = (0.5-difmuk)/s;
+    xx = x*x;
+    fx = -0.5*xx;
+    fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0);
+    if(kflag <= 0) goto S40;
+    goto S60;
+S120:
+/*
+     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
+*/
+//  m = max(1L,(long) (mu));
+    m = (1ULL >= (uint_8)mu) ? 1ULL: (uint_8)mu;
+
+    l = 0;
+    p = exp(-mu);
+    q = p0 = p;
+S130:
+/*
+     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
+*/
+    u = Next();
+    ignpoi = 0;
+    if(u <= p0) return ignpoi;
+/*
+     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
+             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
+             (0.458=PP(9) FOR MU=10)
+*/
+    if(l == 0) goto S150;
+    j = 1;
+//if(u > 0.458) j = min(l,m);
+    if(u > 0.458) j = ((l<=m)? l: m);
+    for(k=j; k<=l; k++) {
+        if(u <= *(pp+k-1)) goto S180;
+    }
+    if(l == 35) goto S130;
+S150:
+/*
+     STEP C. CREATION OF NEW POISSON PROBABILITIES P
+             AND THEIR CUMULATIVES Q=PP(K)
+*/
+    l += 1;
+    for(k=l; k<=35; k++) {
+        p = p*mu/(double)k;
+        q += p;
+        *(pp+k-1) = q;
+        if(u <= q) goto S170;
+    }
+    l = 35;
+    goto S130;
+S170:
+    l = k;
+S180:
+    ignpoi = k;
+    return ignpoi;
+}
+
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+
+r_8 RandomGeneratorInterface::ExpoSimple(void)
+{
+  return -log(1.-Next());
+}
+
+
+static double q_expo[8] = {
+    0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0};
+r_8 RandomGeneratorInterface::ExpoAhrens(void)
+/*
+**********************************************************************
+**********************************************************************
+     (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION                
+**********************************************************************
+**********************************************************************
+                                                                      
+     FOR DETAILS SEE:                                                 
+                                                                      
+               AHRENS, J.H. AND DIETER, U.                            
+               COMPUTER METHODS FOR SAMPLING FROM THE                 
+               EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  
+               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               
+                                                                      
+     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       
+     'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)       
+                                                                      
+     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
+     SUNIF.  The argument IR thus goes away.                          
+                                                                      
+**********************************************************************
+     Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
+     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
+*/
+{
+long i;
+double sexpo,a,u,ustar,umin;
+double *q1 = q_expo;
+    a = 0.0;
+    while((u=Next())==0.);
+    goto S30;
+S20:
+    a += *q1;
+S30:
+    u += u;
+    if(u <= 1.0) goto S20;
+    u -= 1.0;
+    if(u > *q1) goto S60;
+    sexpo = a+u;
+    return sexpo;
+S60:
+    i = 1;
+    ustar = Next();
+    umin = ustar;
+S70:
+    ustar = Next();
+    if(ustar < umin) umin = ustar;
+    i += 1;
+    if(u > *(q_expo+i-1)) goto S70;
+    sexpo = a+umin**q1;
+    return sexpo;
+}
+
+
+}  /* namespace SOPHYA */
+
+
+
+/////////////////////////////////////////////////////////////////
+/*
+**** Remarques sur complex< r_8 > ComplexGaussRan(double sig) ****
+
+--- variables gaussiennes x,y independantes
+x gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2))
+y gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2))
+x,y independants --> pdf f(x,y) = f(x) f(y)
+On a:
+  <x>   = Integrate[x*f(x)]   = Mx
+  <x^2> = Integrate[x^2*f(x)] = Mx^2 + Sx^2
+
+--- On cherche la pdf g(r,t) du module et de la phase
+  x = r cos(t) ,  y = r sin(t)
+  r=sqrt(x^2+y^2 , t=atan2(y,x)
+  (r,t) --> (x,y): le Jacobien = r
+
+  g(r,t) = r f(x,y) = r f(x) f(y)
+         = r/(2Pi Sx Sy) exp(-(x-Mx)^2/(2 Sx^2)) exp(-(y-My)^2/(2 Sy^2))
+
+- Le cas general est complique
+  (cf D.Pelat cours DEA "bruits et signaux" section 4.5)
+
+- Cas ou "Mx = My = 0" et "Sx = Sy = S"
+  c'est la pdf du module et de la phase d'un nombre complexe
+     dont les parties reelles et imaginaires sont independantes
+     et sont distribuees selon des gaussiennes de variance S^2
+  g(r,t) = r/(2Pi S^2) exp(-r^2/(2 S^2))
+  La distribution de "r" est donc:
+    g(r) = Integrate[g(r,t),{t,0,2Pi}]
+         = r/S^2 exp(-r^2/(2 S^2))
+  La distribution de "t" est donc:
+    g(t) = Integrate[g(r,t),{r,0,Infinity}]
+         = 1 / 2Pi  (distribution uniforme sur [0,2Pi[)
+  Les variables aleatoires r,t sont independantes:
+    g(r,t) = g(r) g(t)
+On a:
+  <r>   = Integrate[r*g(r)]   = sqrt(PI/2)*S
+  <r^2> = Integrate[r^2*g(r)] = 2*S^2
+  <r^3> = Integrate[r^3*g(r)] = 3*sqrt(PI/2)*S^3
+  <r^4> = Integrate[r^4*g(r)] = 8*S^4
+
+- Attention:
+La variable complexe "c = x+iy = r*exp(i*t)" ainsi definie verifie:
+              <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 S^2
+Si on veut generer une variable complexe gaussienne telle que
+     <c c*> = s^2 alors il faut prendre S = s/sqrt(2) comme argument
+
+*/
Index: trunk/SophyaLib/BaseTools/randinterf.h
===================================================================
--- trunk/SophyaLib/BaseTools/randinterf.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/randinterf.h	(revision 3602)
@@ -0,0 +1,98 @@
+#ifndef RANDOMGENERATORINTERFACE_H_SEEN
+#define RANDOMGENERATORINTERFACE_H_SEEN
+
+// Classes RandomGeneratorInterface 
+//
+// R. Ansari          (C) UPS+LAL IN2P3/CNRS  2009
+// C. Magneville      (C) DAPNIA/SPP  CEA     2009
+
+
+#include "machdefs.h"
+#include "anydataobj.h"
+#include <complex>
+#include <vector>
+
+namespace SOPHYA {
+
+//! Definition of the interface for Random number generator classes 
+class RandomGeneratorInterface : public AnyDataObj {
+public:
+  RandomGeneratorInterface();
+  virtual ~RandomGeneratorInterface();
+
+  void UseGaussian(uint_2 type=0);
+  void UsePoisson(uint_2 type=0);
+  void UseExpo(uint_2 type=0);
+
+  void GenerateSeedVector(int nseed,vector<uint_2>& seed,int lp=0);
+
+  // --- Le tirage sur une distribution plate
+  /*! \brief Return a random number \b r with a flat distribution  0<=r<1 */ 
+  inline r_8 Flat01() {return Next();}
+  /*! \brief Return a random number \b r with a flat distribution  -1<=r<1 */ 
+  inline r_8 Flatpm1() {return 2.*Next()-1.;}
+
+  // --- Le tirage sur une distribution gaussienne
+  /*! \brief Return a random number following a gaussian distribution (sigma=1, mean=0)*/ 
+  inline r_8 Gaussian() 
+         {
+           if(usegaussian_==0) return GaussianBoxMuller();
+           if(usegaussian_==1) return GaussianSNorm();
+           if(usegaussian_==2) return GaussianPolarBoxMuller();
+           if(usegaussian_==3) return GaussianRatioUnif();
+           if(usegaussian_==4) return GaussianLevaRatioUnif();
+           return GaussianBoxMuller();
+         }
+  virtual r_8 GaussianBoxMuller(); 
+  virtual r_8 GaussianSNorm(); 
+  virtual r_8 GaussianPolarBoxMuller(); 
+  virtual r_8 GaussianRatioUnif(); 
+  virtual r_8 GaussianLevaRatioUnif();
+  /*! \brief Return a random number following a gaussian distribution "sigma", (mean=0)*/ 
+  inline r_8 Gaussian(double sigma) {return sigma*Gaussian();}
+ /*! \brief Return a random number following a gaussian distribution with mean=mu, and sigma */ 
+  inline r_8 Gaussian(double sigma,double mu) {return sigma*Gaussian()+mu;}
+
+  /*! \brief Return a random number following a gaussian tail distribution for x>sdev */ 
+  virtual r_8 GaussianTail(double sdev);
+
+  // --- Le tirage sur une distribution de poisson
+  /*! \brief Return a random number following a poisson distribution with mean mu */ 
+  inline uint_8 Poisson(double mu, double mumax) 
+         {
+           if(usepoisson_==0) return PoissonSimple(mu,mumax);
+           if(usepoisson_==1) return PoissonAhrens(mu);
+           return PoissonSimple(mu,mumax);
+         }
+  virtual uint_8 PoissonSimple(double mu, double mumax=-1);
+  virtual uint_8 PoissonAhrens(double mu);
+
+  // --- Le tirage sur une distribution exponentielle
+  inline r_8 Expo() 
+         {
+           if(useexpo_==0) return ExpoSimple();
+           if(useexpo_==1) return ExpoAhrens();
+           return ExpoSimple();
+         }
+  virtual r_8 ExpoSimple();
+  virtual r_8 ExpoAhrens();
+
+  // --- Le tirage gaussien complexe (cf texte a la fin du .cc)
+  inline complex< r_8 > ComplexGaussRan(void)
+         {return complex< r_8 >(Gaussian(),Gaussian());}
+  inline complex< r_8 > ComplexGaussRan(double sig)
+         {return complex< r_8 >(sig*Gaussian(),sig*Gaussian());}
+  /*! \brief Returns the module of a random complex number generated by ComplexGaussRan */
+  inline double ModComplexGaussRan(double sig=1.)
+         {double r=-log(1.-Next()); return sig*sqrt(2.*r);}
+
+protected:
+  //! Return a random number in [0,1] 
+  virtual r_8 Next();
+  uint_2 usegaussian_, usepoisson_, useexpo_;
+  
+};
+
+} /* namespace SOPHYA */
+
+#endif
Index: trunk/SophyaLib/BaseTools/randr48.cc
===================================================================
--- trunk/SophyaLib/BaseTools/randr48.cc	(revision 3602)
+++ trunk/SophyaLib/BaseTools/randr48.cc	(revision 3602)
@@ -0,0 +1,275 @@
+#include "machdefs.h"
+#include <math.h> 
+#include <stdlib.h>
+#include "thsafeop.h"
+#include "fiondblock.h"
+
+#include "randr48.h"
+
+namespace SOPHYA {
+
+DR48RandGen::DR48RandGen(long int  seed)
+{
+  srand48(seed);
+}
+
+DR48RandGen::~DR48RandGen()
+{
+}
+
+void DR48RandGen::SetSeed(long int  seed)
+{
+  srand48(seed);
+}
+
+void DR48RandGen::SetSeed(uint_2 seed[3])
+{
+  seed48(seed);
+}
+void DR48RandGen::GetSeed(uint_2 seed[3])
+{
+  uint_2 *p, seed_dummy[3] = {0,0,0};
+  p = seed48(seed_dummy);
+  memcpy(seed,p,3*sizeof(uint_2));
+  // on re-initialise a ce qui etait avant
+  seed48(seed);
+}
+
+r_8 DR48RandGen::Next()
+{
+  return drand48();
+}
+
+void DR48RandGen::AutoInit(int lp)
+{
+  vector<uint_2> seed;
+  GenerateSeedVector(0,seed,lp);
+  uint_2 s[3] = {seed[0],seed[1],seed[2]};
+  SetSeed(s);
+}
+
+//------------------------------------------------------------
+//------------------------------------------------------------
+//------------------------------------------------------------
+/*!
+   \class ThSDR48RandGen
+   \ingroup BaseTools
+   \brief Random number generator
+
+   This class is a thread-safe random number generator.
+   Its PPF handler can be used to save the complete state of the class and the underlying
+   random number generator used.
+
+   \sa SOPHYA::ObjFileIO<ThSDR48RandGen>
+
+*/
+
+// Objet statique global pour gestion de lock entre threads
+static ThSafeOp* ths_rand = NULL;
+
+ThSDR48RandGen::ThSDR48RandGen(size_t n, bool tsafe) 
+{
+  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
+  if (tsafe) {   // thread-safe
+    fg_nothrsafe = false;
+    if (n < 1) n = 1024;
+    rseq_.ReSize(n, false);
+    idx_ = n;
+  }
+  else {   // NOT thread-safe
+    fg_nothrsafe = true;
+    idx_ = 1;
+  }
+}
+
+ThSDR48RandGen::ThSDR48RandGen(ThSDR48RandGen const & rg)
+{
+  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
+  if (!rg.fg_nothrsafe) { // thread-safe
+    fg_nothrsafe = false;
+    rseq_.ReSize(rg.rseq_.Size(), false);
+    idx_ = rseq_.Size();
+  }
+  else {   // NOT thread-safe
+    fg_nothrsafe = true;
+    idx_ = 1;
+  }
+}
+
+
+ThSDR48RandGen::~ThSDR48RandGen(void) 
+{
+  // rien a faire 
+}
+
+void ThSDR48RandGen::SetBuffSize(size_t n)
+// redimensionnement du buffer
+{
+  if(fg_nothrsafe) return;
+  if (n < 1) n = 1024;
+  rseq_.ReSize(n, false);
+  idx_ = n;
+}
+
+
+void ThSDR48RandGen::SetSeed(long int seed)
+{
+  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
+  ths_rand->lock();
+  DR48RandGen::SetSeed(seed);
+  ths_rand->unlock();
+  return;
+}
+
+void ThSDR48RandGen::SetSeed(uint_2 seed[3])
+{
+  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
+  ths_rand->lock();
+  SetSeed_P(seed);
+  ths_rand->unlock();
+}
+
+void ThSDR48RandGen::GetSeed(uint_2 seed[3])
+{
+  if (ths_rand == NULL)  ths_rand = new ThSafeOp; 
+  ths_rand->lock();
+  GetSeed_P(seed);
+  ths_rand->unlock();
+  return;
+}
+
+void ThSDR48RandGen::SetSeed_P(uint_2 seed[3])
+{
+  DR48RandGen::SetSeed(seed);
+}
+
+void ThSDR48RandGen::GetSeed_P(uint_2 seed[3])
+{
+  DR48RandGen::GetSeed(seed);
+}
+
+void ThSDR48RandGen::AutoInit(int lp)
+{
+  vector<uint_2> seed;
+  GenerateSeedVector(0,seed,lp);
+  uint_2 s[3] = {seed[0],seed[1],seed[2]};
+  SetSeed(s);
+}
+
+void ThSDR48RandGen::GenSeq(void)
+{
+  ths_rand->lock();
+  for(size_t k=0; k<rseq_.Size(); k++)  rseq_(k) = drand48();
+  ths_rand->unlock();
+  idx_ = 0;
+}
+
+//----------------------------------------------------------
+// Classe pour la gestion de persistance
+// ObjFileIO<DR48RandGen>
+//----------------------------------------------------------
+
+/* --Methode-- */
+DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
+void ObjFileIO<DR48RandGen>::WriteSelf(POutPersist& s) const
+{
+  if (dobj == NULL)
+    throw NullPtrError("ObjFileIO<DR48RandGen>::WriteSelf() dobj=NULL");
+  uint_2 seed[3];
+  dobj->GetSeed(seed);
+  s.Put(seed,3);
+  return;
+}
+
+/* --Methode-- */
+DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
+void ObjFileIO<DR48RandGen>::ReadSelf(PInPersist& s)
+{
+  uint_2 seed[3];
+  s.Get(seed,3);
+  if(dobj == NULL) dobj = new DR48RandGen();
+  dobj->SetSeed(seed);
+  return;
+}
+
+//----------------------------------------------------------
+// Classe pour la gestion de persistance
+// ObjFileIO<ThSDR48RandGen>
+//----------------------------------------------------------
+
+/* --Methode-- */
+DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
+void ObjFileIO<ThSDR48RandGen>::WriteSelf(POutPersist& s) const
+{
+  if (dobj == NULL)
+    throw NullPtrError("ObjFileIO<ThSDR48RandGen>::WriteSelf() dobj=NULL");
+  ths_rand->lock();  // thread-safety
+  uint_4 itab[6];
+  //itab : [0]: version, [1,2,3] = srand48 state/seed  , [4,5] = reserved for future use
+  itab[0] = 1;
+  // On recupere et on ecrit ds le PPF l'etat du generateur aleatoire
+  uint_2 seed_16v[3];
+  dobj->GetSeed_P(seed_16v);
+  for(int i=0; i<3; i++)  itab[i+1] = seed_16v[i];
+  itab[4] = 0;
+  s.Put(itab, 6);
+  uint_8 sz = dobj->rseq_.Size();
+  s.Put(sz);  // Taille du tableau intermediaire
+  uint_8 ix = dobj->idx_;
+  s.Put(ix);  // valeur de l'index 
+
+  if (dobj->rseq_.Size() > 0)  s << dobj->rseq_;  // On ecrit le tableau (NDataBlock) si necessaire
+  ths_rand->unlock(); // thread-safety
+  return;
+}
+
+/* --Methode-- */
+DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
+void ObjFileIO<ThSDR48RandGen>::ReadSelf(PInPersist& s)
+{
+  uint_4 itab[6];
+  //itab : [0]: version, [1,2,3] = srand48 state/seed  , [4] = reserved for future use
+  s.Get(itab, 6);
+  uint_8 sz,ix;
+  s.Get(sz);  // Taille du tableau intermediaire
+  s.Get(ix);  // Taille du tableau intermediaire
+
+  if (dobj == NULL) dobj = new ThSDR48RandGen(sz, (sz>0)?true:false);
+  dobj->idx_ = ix;
+  if (sz > 0) { 
+    s >> dobj->rseq_;  // On lit le tableau (NDataBlock) si necessaire
+    dobj->fg_nothrsafe = false;
+  }
+  else {  // Objet lu est NON thread-safe, taille_tableau rseq_ = 0
+    dobj->fg_nothrsafe = true;
+    if (dobj->rseq_.Size() > 0)  dobj->rseq_.Dealloc(); 
+  }
+  // On initialise l'etat du generateur aleatoire avec les valeurs lues
+  uint_2 seed_16v[3];
+  //NON ? pourquoi faire GetSeed ? : dobj->GetSeed_P(seed_16v);
+  for(int i=0; i<3; i++)  seed_16v[i] = itab[i+1];
+  dobj->SetSeed(seed_16v);
+  return;
+}
+
+// ---------------------------------------------------------
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template ObjFileIO<DR48RandGen>
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class ObjFileIO<DR48RandGen>;
+#endif
+
+// ---------------------------------------------------------
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template ObjFileIO<ThSDR48RandGen>
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class ObjFileIO<ThSDR48RandGen>;
+#endif
+// ---------------------------------------------------------
+
+
+}  /* namespace SOPHYA */
Index: trunk/SophyaLib/BaseTools/randr48.h
===================================================================
--- trunk/SophyaLib/BaseTools/randr48.h	(revision 3602)
+++ trunk/SophyaLib/BaseTools/randr48.h	(revision 3602)
@@ -0,0 +1,98 @@
+#ifndef RANDR48_H_SEEN
+#define RANDR48_H_SEEN
+
+// Classes DR48RandGen
+// Generateur aleatoire utilisant l'aleatoire drand48()
+//
+// R. Ansari          (C) UPS+LAL IN2P3/CNRS  2009
+// C. Magneville      (C) DAPNIA/SPP  CEA     2009
+
+
+#include "machdefs.h"
+#include "objfio.h"
+#include "ndatablock.h"
+#include "randinterf.h"
+
+namespace SOPHYA {
+
+//! Pseudorandom number generator class using dSFMT code 
+class DR48RandGen : public RandomGeneratorInterface {
+
+ public:
+  DR48RandGen(long int seed=12345);
+  virtual ~DR48RandGen();
+
+  virtual void SetSeed(long int seed);
+  virtual void SetSeed(uint_2 seed[3]);
+  virtual void GetSeed(uint_2 seed[3]);
+  virtual void AutoInit(int lp=0);
+
+  friend class ObjFileIO<DR48RandGen> ;
+
+ protected:
+  virtual r_8 Next();
+
+};
+
+//--------------------------------------------------------------------------------
+//--------------------------------------------------------------------------------
+//--------------------------------------------------------------------------------
+
+class ThSDR48RandGen : public DR48RandGen {
+
+ public:
+  ThSDR48RandGen(size_t n=1024,bool tsafe=true);
+  ThSDR48RandGen(ThSDR48RandGen const & rg);
+  virtual ~ThSDR48RandGen();
+
+  void SetBuffSize(size_t n);
+
+  virtual void SetSeed(long int seed);
+  virtual void SetSeed(uint_2 seed[3]);
+  virtual void GetSeed(uint_2 seed[3]);
+  virtual void AutoInit(int lp=0);
+
+  friend class ObjFileIO<ThSDR48RandGen> ;
+
+ protected:
+  // ---- protected data members
+  NDataBlock<r_8>  rseq_;  
+  size_t idx_;
+  bool fg_nothrsafe;  // if true --> NOT thread-safe
+  // ---- protected methods 
+  void GenSeq();
+  inline r_8 Next() 
+    {
+      if (rseq_.Size() == 0)  return drand48(); 
+      else { 
+	if(idx_==rseq_.Size()) GenSeq(); 
+	return(rseq_(idx_++));
+      }
+    }
+  // Non thread-safe version of Init() and GetSeed()
+  void SetSeed_P(uint_2 seed[3]);
+  void GetSeed_P(uint_2 seed[3]);
+
+};
+
+// Classe pour la gestion de persistance PPF :  ObjFileIO<ThSDR48RandGen>
+
+/*! Writes the random generator object state in the POutPersist stream \b os */
+inline POutPersist& operator << (POutPersist& os, ThSDR48RandGen & obj)
+{ ObjFileIO<ThSDR48RandGen> fio(&obj);  fio.Write(os);  return(os); }
+/*! Reads the random generator object state from the PInPersist stream \b is */
+inline PInPersist& operator >> (PInPersist& is, ThSDR48RandGen & obj)
+{ ObjFileIO<ThSDR48RandGen> fio(&obj); is.SkipToNextObject(); fio.Read(is); return(is); }
+
+// Classe pour la gestion de persistance PPF :  ObjFileIO<DR48RandGen>
+
+/*! Writes the random generator object state in the POutPersist stream \b os */
+inline POutPersist& operator << (POutPersist& os, DR48RandGen & obj)
+{ ObjFileIO<DR48RandGen> fio(&obj);  fio.Write(os);  return(os); }
+/*! Reads the random generator object state from the PInPersist stream \b is */
+inline PInPersist& operator >> (PInPersist& is, DR48RandGen & obj)
+{ ObjFileIO<DR48RandGen> fio(&obj); is.SkipToNextObject(); fio.Read(is); return(is); }
+
+} /* namespace SOPHYA */
+
+#endif
Index: trunk/SophyaLib/BaseTools/smakefile
===================================================================
--- trunk/SophyaLib/BaseTools/smakefile	(revision 3601)
+++ trunk/SophyaLib/BaseTools/smakefile	(revision 3602)
@@ -4,8 +4,8 @@
  
 clean:
-	 rm -f $(SOPHYAOBJP)datatype.o $(SOPHYAOBJP)dvlist.o $(SOPHYAOBJP)fiondblock.o $(SOPHYAOBJP)mutyv.o $(SOPHYAOBJP)ndatablock.o $(SOPHYAOBJP)pexceptions.o $(SOPHYAOBJP)ppersist.o $(SOPHYAOBJP)ppfbinstream.o $(SOPHYAOBJP)ppftpointerio.o $(SOPHYAOBJP)rawstream.o $(SOPHYAOBJP)sophyainit.o $(SOPHYAOBJP)stsrand.o $(SOPHYAOBJP)timestamp.o $(SOPHYAOBJP)gnumd5.o $(SOPHYAOBJP)srandgen.o
+	 rm -f $(SOPHYAOBJP)datatype.o $(SOPHYAOBJP)dvlist.o $(SOPHYAOBJP)fiondblock.o $(SOPHYAOBJP)mutyv.o $(SOPHYAOBJP)ndatablock.o $(SOPHYAOBJP)pexceptions.o $(SOPHYAOBJP)ppersist.o $(SOPHYAOBJP)ppfbinstream.o $(SOPHYAOBJP)ppftpointerio.o $(SOPHYAOBJP)randfmt.o $(SOPHYAOBJP)randinterf.o $(SOPHYAOBJP)randr48.o $(SOPHYAOBJP)rawstream.o $(SOPHYAOBJP)sophyainit.o $(SOPHYAOBJP)stsrand.o $(SOPHYAOBJP)timestamp.o $(SOPHYAOBJP)dSFMT.o $(SOPHYAOBJP)gnumd5.o $(SOPHYAOBJP)srandgen.o
 	rm -f $(SOPHYALIBP)libBaseTools.a
  
-$(SOPHYALIBP)libBaseTools.a : $(SOPHYAOBJP)datatype.o $(SOPHYAOBJP)dvlist.o $(SOPHYAOBJP)fiondblock.o $(SOPHYAOBJP)gnumd5.o $(SOPHYAOBJP)mutyv.o $(SOPHYAOBJP)ndatablock.o $(SOPHYAOBJP)pexceptions.o $(SOPHYAOBJP)ppersist.o $(SOPHYAOBJP)ppfbinstream.o $(SOPHYAOBJP)ppftpointerio.o $(SOPHYAOBJP)rawstream.o $(SOPHYAOBJP)sophyainit.o $(SOPHYAOBJP)srandgen.o $(SOPHYAOBJP)stsrand.o $(SOPHYAOBJP)timestamp.o
+$(SOPHYALIBP)libBaseTools.a : $(SOPHYAOBJP)datatype.o $(SOPHYAOBJP)dSFMT.o $(SOPHYAOBJP)dvlist.o $(SOPHYAOBJP)fiondblock.o $(SOPHYAOBJP)gnumd5.o $(SOPHYAOBJP)mutyv.o $(SOPHYAOBJP)ndatablock.o $(SOPHYAOBJP)pexceptions.o $(SOPHYAOBJP)ppersist.o $(SOPHYAOBJP)ppfbinstream.o $(SOPHYAOBJP)ppftpointerio.o $(SOPHYAOBJP)randfmt.o $(SOPHYAOBJP)randinterf.o $(SOPHYAOBJP)randr48.o $(SOPHYAOBJP)rawstream.o $(SOPHYAOBJP)sophyainit.o $(SOPHYAOBJP)srandgen.o $(SOPHYAOBJP)stsrand.o $(SOPHYAOBJP)timestamp.o
 	$(AR) $(ARFLAGS) $@ $? 
 	touch $(SOPHYAINCP)/SophyaConfInfo/libsophya.objlist
@@ -59,4 +59,21 @@
 	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  ppftpointerio.cc
  
+$(SOPHYAOBJP)randfmt.o: randfmt.cc randfmt.h \
+  $(SOPHYAINCP)machdefs.h objfio.h anydataobj.h \
+  ppersist.h pexceptions.h gnumd5.h ppfbinstream.h rawstream.h \
+  ppfnametag.h dsfmtflags.h dSFMT.h randinterf.h
+	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  randfmt.cc
+ 
+$(SOPHYAOBJP)randinterf.o: randinterf.cc sopnamsp.h \
+  $(SOPHYAINCP)machdefs.h pexceptions.h \
+  randinterf.h anydataobj.h
+	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  randinterf.cc
+ 
+$(SOPHYAOBJP)randr48.o: randr48.cc $(SOPHYAINCP)machdefs.h \
+  thsafeop.h pexceptions.h fiondblock.h ppersist.h gnumd5.h \
+  ppfbinstream.h rawstream.h ppfnametag.h ndatablock.h anydataobj.h \
+  randr48.h objfio.h randinterf.h
+	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  randr48.cc
+ 
 $(SOPHYAOBJP)rawstream.o: rawstream.cc sopnamsp.h rawstream.h \
   $(SOPHYAINCP)machdefs.h pexceptions.h
@@ -67,5 +84,6 @@
   pexceptions.h ppersist.h gnumd5.h ppfbinstream.h rawstream.h \
   ppfnametag.h fiondblock.h ndatablock.h anydataobj.h dvlist.h objfio.h \
-  mutyv.h timestamp.h stsrand.h fiosegdb.h segdatablock.h thsafeop.h \
+  mutyv.h timestamp.h stsrand.h randr48.h randinterf.h randfmt.h \
+  dsfmtflags.h dSFMT.h fiosegdb.h segdatablock.h thsafeop.h \
   ppftpointerio.h ppfwrapstlv.h sversion.h
 	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  sophyainit.cc
@@ -82,4 +100,7 @@
 	$(CXXCOMPILE) $(CXXTEMPFLG) -o $@  timestamp.cc
  
+$(SOPHYAOBJP)dSFMT.o: dSFMT.c dsfmtflags.h dSFMT-params.h dSFMT.h dSFMT-params19937.h
+	$(CCOMPILE) -o $@  dSFMT.c
+ 
 $(SOPHYAOBJP)gnumd5.o: gnumd5.c $(SOPHYAINCP)machdefs.h \
   bithelp.h gnumd5.h
Index: trunk/SophyaLib/BaseTools/sophyainit.cc
===================================================================
--- trunk/SophyaLib/BaseTools/sophyainit.cc	(revision 3601)
+++ trunk/SophyaLib/BaseTools/sophyainit.cc	(revision 3602)
@@ -16,4 +16,6 @@
 
 #include "stsrand.h"
+#include "randr48.h"
+#include "randfmt.h"
 
 #include "fiosegdb.h" 
@@ -101,4 +103,10 @@
   PPRegister(ObjFileIO<RandomGenerator>);
   DObjRegister(ObjFileIO<RandomGenerator>, RandomGenerator);
+  PPRegister(ObjFileIO<DR48RandGen>);
+  DObjRegister(ObjFileIO<DR48RandGen>,DR48RandGen );
+  PPRegister(ObjFileIO<ThSDR48RandGen>);
+  DObjRegister(ObjFileIO<ThSDR48RandGen>, ThSDR48RandGen);
+  PPRegister(ObjFileIO<FMTRandGen>);
+  DObjRegister(ObjFileIO<FMTRandGen>, FMTRandGen);
 
   // Enregistrement des handlers PPF pour les SegDataBlock<T> 
