Changeset 4019 in Sophya for trunk/SophyaLib/BaseTools/randinterf.cc

Timestamp:

Sep 24, 2011, 6:31:08 PM (14 years ago)

Author:

cmv

Message:

ajout generateurs tinymt32+64, cmv 24/09/2011

File:

• trunk/SophyaLib/BaseTools/randinterf.cc (modified) (9 diffs)

trunk/SophyaLib/BaseTools/randinterf.cc

@@ -88 +88 @@
 /////////////////////////////////////////////////////////////////////////
 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
+// renvoie un vecteur de 3+nseed entiers 16 bits
+// [0 - 2] = codage sur 3*16 = 48 bits du nombre (melange) de microsec depuis l'origine
 // [3 -> 3+ngene-1] = entiers aleatoires (poor man generator)
 //
@@ -112 +112 @@
   // 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;
+  if(lp>1) cout<<".since orig: "<<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;
@@ -125 +125 @@
   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]);
+    else if(ip%3==2) swap(b[ip],b[16+ip]);
   }
   if(lp>2) {
@@ -139 +139 @@
     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;
+  if(lp>0) cout<<"seed_time[0-2]: "<<seed[0]<<" "<<seed[1]<<" "<<seed[2]<<endl;
+  // On recree tmicro70 avec les bits melanges
+  tmicro70 = uint_8(seed[0]) | (uint_8(seed[1])<<16) | (uint_8(seed[2])<<32);
 
   // ---
@@ -160 +162 @@
   // the requirements. The constants below are arbitrary otherwise
   uint_4 seed0 = uint_4(tmicro70&0xFFFFFFFFULL);
-  if(lp>2) cout<<"seed0(time): "<<seed0<<endl;
+  if(lp>2) cout<<"seed0_time_init_t88: "<<seed0<<endl;
   uint_4 state_s1, state_s2, state_s3;
   state_s1 = 1243598713U ^ seed0; if (state_s1 <  2) state_s1 = 1243598713U;
@@ -177 +179 @@
     uint_2 s = uint_2( (s1^s2^s3)&0xFFFFU );
     seed.push_back(s);
-    if(lp>0) cout<<"seed(t88): "<<seed[3+nfill]<<endl;
+    if(lp>0) cout<<"seed_t88["<<nfill+3<<"]: "<<seed[3+nfill]<<endl;
     nfill++;
   }
@@ -356 +358 @@
 r_8 RandomGeneratorInterface::GaussianPolarBoxMuller()
 {
-double x1,x2,w;
-do {
+ 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);
+ } while ( w >= 1.0 || w==0. );
+ return x1 * sqrt(-2.0*log(w)/w);
 }
 
@@ -368 +370 @@
 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;
+ 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;
+ double u,v,x,y,q;
  do {
    u = 1.-Next();  // in ]0,1]
@@ -525 +527 @@
    * This implementation does one-sided upper-tailed deviates.
    */
-
-  if (s < slim)
-    {
-      /* 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;
-    }
+  if(s < slim) {
+    /* 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;
+  }
 }