src/engine/Plugin/aleas.cxx

   1 // --- C++ ---
   2 // --- coding: latin_1 ---
   3 //
   4 //    File
   5 //      creation : 2007-02-21.09.50.46
   6 //      revision : $Id$
   7 //
   8 //    Copyright © 2007 Commissariat à l'Energie Atomique
   9 //      par Gilles ARNAUD (DM2S/SFME/LETR)
  10 //        C.E. Saclay; Bat 454; 91191 GIF/YVETTE CEDEX; France
  11 //        Tel: 01 69 08 38 86; Fax : 33 1 69 08 85 68
  12 //        Gilles.Arnaud@cea.fr
  13 //
  14 //    Object
  15 //      tirage aleatoire
  16 //
  17 //___________________________________________________________________
  18
  19 #include <time.h>
  20 #include <sys/types.h>
  21 #include <unistd.h>
  22 #include <math.h>
  23
  24 // #include "mt19937ar.h
  25 extern void init_by_array(unsigned long [], int);
  26 extern double genrand_real1(void);
  27 extern double genrand_real3(void);
  28
  29 #include "aleas.hxx"
  30
  31 #define NOALEAS
  32
  33 static void initrand(void)
  34 {
  35      static long done=0;
  36      unsigned long vec[]={
  37              31081996, 21012006, 17921963, 0,
  38              11101998, 2112003, 25111964, 0
  39          };
  40
  41      if (! done) {
  42 #ifdef NOALEAS
  43          vec[3] = vec[7] = 2082007;
  44 #else
  45          vec[3] = (unsigned long) getpid();
  46          vec[7] = (unsigned long) time(NULL);
  47 #endif
  48          init_by_array(vec, 8);
  49          done += 1;
  50      }
  51 }
  52 // utilitaire
  53 /*static void initrand(void)
  54 {
  55     static long done=0;
  56     unsigned long vec[]={
  57             31081996, 21012006, 17921963, 0,
  58             11101998, 2112003, 25111964, 0
  59         };
  60
  61     if (! done) {
  62         vec[3] = (unsigned long) getpid();
  63         vec[7] = (unsigned long) time(NULL);
  64         init_by_array(vec, 8);
  65         done += 1;
  66     }
  67 }*/
  68
  69 static double randGauss(void)
  70 {
  71     static double   next;
  72     static int      flag=0;
  73     double          v2, d, fac;
  74
  75     if (flag) {
  76         flag = 0;
  77         return next;
  78     }
  79     else {
  80         do {
  81             next = 2.0*genrand_real3() - 1.0;
  82             v2 = 2.0*genrand_real3() - 1.0;
  83             d = next*next + v2*v2;
  84         } while (d >= 1.0 || d == 0.0);
  85         fac = sqrt(-2.0*log(d)/d);
  86         next *= fac;
  87         flag = 1;
  88         return v2 * fac;
  89     }
  90
  91 }
  92
  93 // class Aleatoire
  94 Aleatoire::Aleatoire(long sz)
  95 {
  96     size = sz;
  97     initrand();
  98 }
  99
 100 void Aleatoire::fill(std::vector<double> &ret)
 101 {
 102     int     i;
 103
 104     for (i=0; i<size; i++) ret[i] = tire();
 105 }
 106
 107 std::vector<double> *Aleatoire::gen()
 108 {
 109     std::vector<double> *ret;
 110
 111     ret = new std::vector<double>(size);
 112     fill(*ret);
 113     return ret;
 114 }
 115
 116 // class Cube
 117 double Cube::tire()
 118 {
 119     return genrand_real1();
 120 }
 121
 122 // class NormalPositif
 123 double NormalePositif::tire()
 124 {
 125     return randGauss() * 0.25 + 0.5 ;
 126 }
 127
 128 // class Normal
 129 double Normale::tire()
 130 {
 131     return randGauss();
 132 }
 133
 134 // class Sphere
 135 void Sphere::fill(std::vector<double> &ret)
 136 {
 137     long    i;
 138     double  cum, r;
 139
 140     Normale::fill(ret);
 141     for (cum=0, i=0; i<size; i++)
 142         cum += ret[i] * ret[i];
 143     cum = sqrt(cum);
 144     r = pow(genrand_real1(), size);
 145     for (i=0; i<size; i++)
 146         ret[i] *= cum * r;
 147 }
 148
 149 // class SpherePositif
 150 void SpherePositif::fill(std::vector<double> &ret)
 151 {
 152     long    i;
 153
 154     Sphere::fill(ret);
 155     for (i=0; i<size; i++)
 156         ret[i] = fabs(ret[i]);
 157 }
 158