libgfortran/intrinsics/random.c

   1 /* Implementation of the RANDOM intrinsics
   2    Copyright 2002, 2004 Free Software Foundation, Inc.
   3    Contributed by Lars Segerlund <seger@linuxmail.org>
   4
   5   The algorithm was taken from the paper :
   6
   7         Mersenne Twister:       623-dimensionally equidistributed
   8                                 uniform pseudorandom generator.
   9
  10         by:     Makoto Matsumoto
  11                 Takuji Nishimura
  12
  13         Which appeared in the: ACM Transactions on Modelling and Computer
  14         Simulations: Special Issue on Uniform Random Number
  15         Generation. ( Early in 1998 ).
  16
  17 This file is part of the GNU Fortran 95 runtime library (libgfortran).
  18
  19 Libgfortran is free software; you can redistribute it and/or
  20 modify it under the terms of the GNU Lesser General Public
  21 License as published by the Free Software Foundation; either
  22 version 2.1 of the License, or (at your option) any later version.
  23
  24 Ligbfortran is distributed in the hope that it will be useful,
  25 but WITHOUT ANY WARRANTY; without even the implied warranty of
  26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  27 GNU Lesser General Public License for more details.
  28
  29 You should have received a copy of the GNU Lesser General Public
  30 License along with libgfor; see the file COPYING.LIB.  If not,
  31 write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  32 Boston, MA 02111-1307, USA.  */
  33
  34 #include "config.h"
  35 #include <stdio.h>
  36 #include <stdlib.h>
  37 #include <sys/types.h>
  38 #include <sys/stat.h>
  39 #include <fcntl.h>
  40
  41 #ifdef HAVE_UNISTD_H
  42 #include <unistd.h>
  43 #endif
  44
  45 #include "libgfortran.h"
  46
  47 /*Use the 'big' generator by default ( period -> 2**19937 ).  */
  48
  49 #define MT19937
  50
  51 /* Define the necessary constants for the algorithm.  */
  52
  53 #ifdef  MT19937
  54 enum constants
  55 {
  56   N = 624, M = 397, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
  57 };
  58 #define M_A     0x9908B0DF
  59 #define T_B     0x9D2C5680
  60 #define T_C     0xEFC60000
  61 #else
  62 enum constants
  63 {
  64   N = 351, M = 175, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
  65 };
  66 #define M_A     0xE4BD75F5
  67 #define T_B     0x655E5280
  68 #define T_C     0xFFD58000
  69 #endif
  70
  71 static int i = N;
  72 static unsigned int seed[N];
  73
  74 /* This is the routine which handles the seeding of the generator,
  75    and also reading and writing of the seed.  */
  76
  77 void
  78 random_seed (GFC_INTEGER_4 * size, const gfc_array_i4 * put,
  79              const gfc_array_i4 * get)
  80 {
  81   /* Initialize the seed in system dependent manner.  */
  82   if (get == NULL && put == NULL && size == NULL)
  83     {
  84       int fd;
  85       fd = open ("/dev/urandom", O_RDONLY);
  86       if (fd == 0)
  87         {
  88           /* We dont have urandom.  */
  89           GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed;
  90           for (i = 0; i < N; i++)
  91             {
  92               s = s * 29943829 - 1;
  93               seed[i] = s;
  94             }
  95         }
  96       else
  97         {
  98           /* Using urandom, might have a length issue.  */
  99           read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N);
 100           close (fd);
 101         }
 102       return;
 103     }
 104
 105   /* Return the size of the seed */
 106   if (size != NULL)
 107     {
 108       *size = N;
 109       return;
 110     }
 111
 112   /* if we have gotten to this pount we have a get or put
 113    * now we check it the array fulfills the demands in the standard .
 114    */
 115
 116   /* Set the seed to PUT data */
 117   if (put != NULL)
 118     {
 119       /* if the rank of the array is not 1 abort */
 120       if (GFC_DESCRIPTOR_RANK (put) != 1)
 121         abort ();
 122
 123       /* if the array is too small abort */
 124       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N)
 125         abort ();
 126
 127       /* If this is the case the array is a temporary */
 128       if (put->dim[0].stride == 0)
 129         return;
 130
 131       /*  This code now should do correct strides. */
 132       for (i = 0; i < N; i++)
 133         seed[i] = put->data[i * put->dim[0].stride];
 134     }
 135
 136   /* Return the seed to GET data */
 137   if (get != NULL)
 138     {
 139       /* if the rank of the array is not 1 abort */
 140       if (GFC_DESCRIPTOR_RANK (get) != 1)
 141         abort ();
 142
 143       /* if the array is too small abort */
 144       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N)
 145         abort ();
 146
 147       /* If this is the case the array is a temporary */
 148       if (get->dim[0].stride == 0)
 149         return;
 150
 151       /*  This code now should do correct strides. */
 152       for (i = 0; i < N; i++)
 153         get->data[i * get->dim[0].stride] = seed[i];
 154     }
 155 }
 156
 157 /* Here is the internal routine which generates the random numbers
 158    in 'batches' based upon the need for a new batch.
 159    It's an integer based routine known as 'Mersenne Twister'.
 160    This implementation still lacks 'tempering' and a good verification,
 161    but gives very good metrics.  */
 162
 163 static void
 164 random_generate (void)
 165 {
 166   /* 32 bits.  */
 167   GFC_UINTEGER_4 y;
 168
 169   /* Generate batch of N.  */
 170   int k, m;
 171   for (k = 0, m = M; k < N - 1; k++)
 172     {
 173       y = (seed[k] & (-1 << R)) | (seed[k + 1] & ((1u << R) - 1));
 174       seed[k] = seed[m] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
 175       if (++m >= N)
 176         m = 0;
 177     }
 178
 179   y = (seed[N - 1] & (-1 << R)) | (seed[0] & ((1u << R) - 1));
 180   seed[N - 1] = seed[M - 1] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
 181   i = 0;
 182 }
 183
 184 /* A routine to return a REAL(KIND=4).  */
 185
 186 #define random_r4 prefix(random_r4)
 187 void
 188 random_r4 (GFC_REAL_4 * harv)
 189 {
 190   /* Regenerate if we need to.  */
 191   if (i >= N)
 192     random_generate ();
 193
 194   /* Convert uint32 to REAL(KIND=4).  */
 195   *harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
 196                         (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
 197 }
 198
 199 /* A routine to return a REAL(KIND=8).  */
 200
 201 #define random_r8 prefix(random_r8)
 202 void
 203 random_r8 (GFC_REAL_8 * harv)
 204 {
 205   /* Regenerate if we need to, may waste one 32-bit value.  */
 206   if ((i + 1) >= N)
 207     random_generate ();
 208
 209   /* Convert two uint32 to a REAL(KIND=8).  */
 210   *harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
 211           (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
 212   i += 2;
 213 }
 214
 215 /* Code to handle arrays will follow here.  */
 216
 217 /* REAL(KIND=4) REAL array.  */
 218
 219 #define arandom_r4 prefix(arandom_r4)
 220 void
 221 arandom_r4 (gfc_array_r4 * harv)
 222 {
 223   index_type count[GFC_MAX_DIMENSIONS - 1];
 224   index_type extent[GFC_MAX_DIMENSIONS - 1];
 225   index_type stride[GFC_MAX_DIMENSIONS - 1];
 226   index_type stride0;
 227   index_type dim;
 228   GFC_REAL_4 *dest;
 229   int n;
 230
 231   dest = harv->data;
 232
 233   if (harv->dim[0].stride == 0)
 234     harv->dim[0].stride = 1;
 235
 236   dim = GFC_DESCRIPTOR_RANK (harv);
 237
 238   for (n = 0; n < dim; n++)
 239     {
 240       count[n] = 0;
 241       stride[n] = harv->dim[n].stride;
 242       extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
 243       if (extent[n] <= 0)
 244         return;
 245     }
 246
 247   stride0 = stride[0];
 248
 249   while (dest)
 250     {
 251       /* Set the elements.  */
 252
 253       /* regenerate if we need to */
 254       if (i >= N)
 255         random_generate ();
 256
 257       /* Convert uint32 to float in a hopefully g95 compiant manner */
 258       *dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
 259                             (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
 260
 261       /* Advance to the next element.  */
 262       dest += stride0;
 263       count[0]++;
 264       /* Advance to the next source element.  */
 265       n = 0;
 266       while (count[n] == extent[n])
 267         {
 268           /* When we get to the end of a dimension,
 269              reset it and increment
 270              the next dimension.  */
 271           count[n] = 0;
 272           /* We could precalculate these products,
 273              but this is a less
 274              frequently used path so proabably not worth it.  */
 275           dest -= stride[n] * extent[n];
 276           n++;
 277           if (n == dim)
 278             {
 279               dest = NULL;
 280               break;
 281             }
 282           else
 283             {
 284               count[n]++;
 285               dest += stride[n];
 286             }
 287         }
 288     }
 289 }
 290
 291 /* REAL(KIND=8) array.  */
 292
 293 #define arandom_r8 prefix(arandom_r8)
 294 void
 295 arandom_r8 (gfc_array_r8 * harv)
 296 {
 297   index_type count[GFC_MAX_DIMENSIONS - 1];
 298   index_type extent[GFC_MAX_DIMENSIONS - 1];
 299   index_type stride[GFC_MAX_DIMENSIONS - 1];
 300   index_type stride0;
 301   index_type dim;
 302   GFC_REAL_8 *dest;
 303   int n;
 304
 305   dest = harv->data;
 306
 307   if (harv->dim[0].stride == 0)
 308     harv->dim[0].stride = 1;
 309
 310   dim = GFC_DESCRIPTOR_RANK (harv);
 311
 312   for (n = 0; n < dim; n++)
 313     {
 314       count[n] = 0;
 315       stride[n] = harv->dim[n].stride;
 316       extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
 317       if (extent[n] <= 0)
 318         return;
 319     }
 320
 321   stride0 = stride[0];
 322
 323   while (dest)
 324     {
 325       /* Set the elements.  */
 326
 327       /* regenerate if we need to, may waste one 32-bit value */
 328       if ((i + 1) >= N)
 329         random_generate ();
 330
 331       /* Convert two uint32 to a REAL(KIND=8).  */
 332       *dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
 333               (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
 334       i += 2;
 335
 336       /* Advance to the next element.  */
 337       dest += stride0;
 338       count[0]++;
 339       /* Advance to the next source element.  */
 340       n = 0;
 341       while (count[n] == extent[n])
 342         {
 343           /* When we get to the end of a dimension,
 344              reset it and increment
 345              the next dimension.  */
 346           count[n] = 0;
 347           /* We could precalculate these products,
 348              but this is a less
 349              frequently used path so proabably not worth it.  */
 350           dest -= stride[n] * extent[n];
 351           n++;
 352           if (n == dim)
 353             {
 354               dest = NULL;
 355               break;
 356             }
 357           else
 358             {
 359               count[n]++;
 360               dest += stride[n];
 361             }
 362         }
 363     }
 364 }
 365