src/engine/Plugin/aleas.cxx

   1 // Copyright (C) 2006-2014  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include <time.h>
  21 #include <sys/types.h>
  22 #ifdef WIN32
  23 #else
  24 #include <unistd.h>
  25 #endif
  26 #include <math.h>
  27
  28 // #include "mt19937ar.h
  29 extern void init_by_array(unsigned long [], int);
  30 extern double genrand_real1(void);
  31 extern double genrand_real3(void);
  32
  33 #include "aleas.hxx"
  34
  35 #define NOALEAS
  36
  37 static void initrand(void)
  38 {
  39      static long done=0;
  40      unsigned long vec[]={
  41              31081996, 21012006, 17921963, 0,
  42              11101998, 2112003, 25111964, 0
  43          };
  44
  45      if (! done) {
  46 #ifdef NOALEAS
  47          vec[3] = vec[7] = 2082007;
  48 #else
  49          vec[3] = (unsigned long) getpid();
  50          vec[7] = (unsigned long) time(NULL);
  51 #endif
  52          init_by_array(vec, 8);
  53          done += 1;
  54      }
  55 }
  56 // utilitaire
  57 /*static void initrand(void)
  58 {
  59     static long done=0;
  60     unsigned long vec[]={
  61             31081996, 21012006, 17921963, 0,
  62             11101998, 2112003, 25111964, 0
  63         };
  64
  65     if (! done) {
  66         vec[3] = (unsigned long) getpid();
  67         vec[7] = (unsigned long) time(NULL);
  68         init_by_array(vec, 8);
  69         done += 1;
  70     }
  71 }*/
  72
  73 static double randGauss(void)
  74 {
  75     static double   next;
  76     static int      flag=0;
  77     double          v2, d, fac;
  78
  79     if (flag) {
  80         flag = 0;
  81         return next;
  82     }
  83     else {
  84         do {
  85             next = 2.0*genrand_real3() - 1.0;
  86             v2 = 2.0*genrand_real3() - 1.0;
  87             d = next*next + v2*v2;
  88         } while (d >= 1.0 || d == 0.0);
  89         fac = sqrt(-2.0*log(d)/d);
  90         next *= fac;
  91         flag = 1;
  92         return v2 * fac;
  93     }
  94
  95 }
  96
  97 // class Aleatoire
  98 Aleatoire::Aleatoire(long sz)
  99 {
 100     size = sz;
 101     initrand();
 102 }
 103
 104 void Aleatoire::fill(std::vector<double> &ret)
 105 {
 106     int     i;
 107
 108     for (i=0; i<size; i++) ret[i] = tire();
 109 }
 110
 111 std::vector<double> *Aleatoire::gen()
 112 {
 113     std::vector<double> *ret;
 114
 115     ret = new std::vector<double>(size);
 116     fill(*ret);
 117     return ret;
 118 }
 119
 120 // class Cube
 121 double Cube::tire()
 122 {
 123     return genrand_real1();
 124 }
 125
 126 // class NormalPositif
 127 double NormalePositif::tire()
 128 {
 129     return randGauss() * 0.25 + 0.5 ;
 130 }
 131
 132 // class Normal
 133 double Normale::tire()
 134 {
 135     return randGauss();
 136 }
 137
 138 // class Sphere
 139 void Sphere::fill(std::vector<double> &ret)
 140 {
 141     long    i;
 142     double  cum, r;
 143
 144     Normale::fill(ret);
 145     for (cum=0, i=0; i<size; i++)
 146         cum += ret[i] * ret[i];
 147     cum = sqrt(cum);
 148     r = pow(genrand_real1(), size);
 149     for (i=0; i<size; i++)
 150         ret[i] *= cum * r;
 151 }
 152
 153 // class SpherePositif
 154 void SpherePositif::fill(std::vector<double> &ret)
 155 {
 156     long    i;
 157
 158     Sphere::fill(ret);
 159     for (i=0; i<size; i++)
 160         ret[i] = fabs(ret[i]);
 161 }
 162