libgfortran/intrinsics/random.c

   1 /* Implementation of the RANDOM intrinsics
   2    Copyright 2002 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 #include <assert.h>
  41 #include "libgfortran.h"
  42
  43 /*Use the 'big' generator by default ( period -> 2**19937 ).  */
  44
  45 #define MT19937
  46
  47 /* Define the necessary constants for the algorithm.  */
  48
  49 #ifdef  MT19937
  50 enum constants
  51 {
  52   N = 624, M = 397, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
  53 };
  54 #define M_A     0x9908B0DF
  55 #define T_B     0x9D2C5680
  56 #define T_C     0xEFC60000
  57 #else
  58 enum constants
  59 {
  60   N = 351, M = 175, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
  61 };
  62 #define M_A     0xE4BD75F5
  63 #define T_B     0x655E5280
  64 #define T_C     0xFFD58000
  65 #endif
  66
  67 static int i = N;
  68 static unsigned int seed[N];
  69
  70 /* This is the routine which handles the seeding of the generator,
  71    and also reading and writing of the seed.  */
  72
  73 #define random_seed prefix(random_seed)
  74 void
  75 random_seed (GFC_INTEGER_4 * size, const gfc_array_i4 * put,
  76              const gfc_array_i4 * get)
  77 {
  78   /* Initialize the seed in system dependent manner.  */
  79   if (get == NULL && put == NULL && size == NULL)
  80     {
  81       int fd;
  82       fd = open ("/dev/urandom", O_RDONLY);
  83       if (fd == 0)
  84         {
  85           /* We dont have urandom.  */
  86           GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed;
  87           for (i = 0; i < N; i++)
  88             {
  89               s = s * 29943829 - 1;
  90               seed[i] = s;
  91             }
  92         }
  93       else
  94         {
  95           /* Using urandom, might have a length issue.  */
  96           read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N);
  97           close (fd);
  98         }
  99       return;
 100     }
 101
 102   /* Return the size of the seed */
 103   if (size != NULL)
 104     {
 105       *size = N;
 106       return;
 107     }
 108
 109   /* if we have gotten to this pount we have a get or put
 110    * now we check it the array fulfills the demands in the standard .
 111    */
 112
 113   /* Set the seed to PUT data */
 114   if (put != NULL)
 115     {
 116       /* if the rank of the array is not 1 abort */
 117       if (GFC_DESCRIPTOR_RANK (put) != 1)
 118         abort ();
 119
 120       /* if the array is too small abort */
 121       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N)
 122         abort ();
 123
 124       /* If this is the case the array is a temporary */
 125       if (get->dim[0].stride == 0)
 126         return;
 127
 128       /*  This code now should do correct strides. */
 129       for (i = 0; i < N; i++)
 130         seed[i] = put->data[i * put->dim[0].stride];
 131     }
 132
 133   /* Return the seed to GET data */
 134   if (get != NULL)
 135     {
 136       /* if the rank of the array is not 1 abort */
 137       if (GFC_DESCRIPTOR_RANK (get) != 1)
 138         abort ();
 139
 140       /* if the array is too small abort */
 141       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N)
 142         abort ();
 143
 144       /* If this is the case the array is a temporary */
 145       if (get->dim[0].stride == 0)
 146         return;
 147
 148       /*  This code now should do correct strides. */
 149       for (i = 0; i < N; i++)
 150         get->data[i * get->dim[0].stride] = seed[i];
 151     }
 152 }
 153
 154 /* Here is the internal routine which generates the random numbers
 155    in 'batches' based upon the need for a new batch.
 156    It's an integer based routine known as 'Mersenne Twister'.
 157    This implementation still lacks 'tempering' and a good verification,
 158    but gives very good metrics.  */
 159
 160 static void
 161 random_generate (void)
 162 {
 163   /* 32 bits.  */
 164   GFC_UINTEGER_4 y;
 165
 166   /* Generate batch of N.  */
 167   int k, m;
 168   for (k = 0, m = M; k < N - 1; k++)
 169     {
 170       y = (seed[k] & (-1 << R)) | (seed[k + 1] & ((1u << R) - 1));
 171       seed[k] = seed[m] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
 172       if (++m >= N)
 173         m = 0;
 174     }
 175
 176   y = (seed[N - 1] & (-1 << R)) | (seed[0] & ((1u << R) - 1));
 177   seed[N - 1] = seed[M - 1] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
 178   i = 0;
 179 }
 180
 181 /* A routine to return a REAL(KIND=4).  */
 182
 183 #define random_r4 prefix(random_r4)
 184 void
 185 random_r4 (GFC_REAL_4 * harv)
 186 {
 187   /* Regenerate if we need to.  */
 188   if (i >= N)
 189     random_generate ();
 190
 191   /* Convert uint32 to REAL(KIND=4).  */
 192   *harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
 193                         (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
 194 }
 195
 196 /* A routine to return a REAL(KIND=8).  */
 197
 198 #define random_r8 prefix(random_r8)
 199 void
 200 random_r8 (GFC_REAL_8 * harv)
 201 {
 202   /* Regenerate if we need to, may waste one 32-bit value.  */
 203   if ((i + 1) >= N)
 204     random_generate ();
 205
 206   /* Convert two uint32 to a REAL(KIND=8).  */
 207   *harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
 208           (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
 209   i += 2;
 210 }
 211
 212 /* Code to handle arrays will follow here.  */
 213
 214 /* REAL(KIND=4) REAL array.  */
 215
 216 #define arandom_r4 prefix(arandom_r4)
 217 void
 218 arandom_r4 (gfc_array_r4 * harv)
 219 {
 220   index_type count[GFC_MAX_DIMENSIONS - 1];
 221   index_type extent[GFC_MAX_DIMENSIONS - 1];
 222   index_type stride[GFC_MAX_DIMENSIONS - 1];
 223   index_type stride0;
 224   index_type dim;
 225   GFC_REAL_4 *dest;
 226   int n;
 227
 228   dest = harv->data;
 229
 230   if (harv->dim[0].stride == 0)
 231     harv->dim[0].stride = 1;
 232
 233   dim = GFC_DESCRIPTOR_RANK (harv);
 234
 235   for (n = 0; n < dim; n++)
 236     {
 237       count[n] = 0;
 238       stride[n] = harv->dim[n].stride;
 239       extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
 240       if (extent[n] <= 0)
 241         return;
 242     }
 243
 244   stride0 = stride[0];
 245
 246   while (dest)
 247     {
 248       /* Set the elements.  */
 249
 250       /* regenerate if we need to */
 251       if (i >= N)
 252         random_generate ();
 253
 254       /* Convert uint32 to float in a hopefully g95 compiant manner */
 255       *dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
 256                             (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
 257
 258       /* Advance to the next element.  */
 259       dest += stride0;
 260       count[0]++;
 261       /* Advance to the next source element.  */
 262       n = 0;
 263       while (count[n] == extent[n])
 264         {
 265           /* When we get to the end of a dimension,
 266              reset it and increment
 267              the next dimension.  */
 268           count[n] = 0;
 269           /* We could precalculate these products,
 270              but this is a less
 271              frequently used path so proabably not worth it.  */
 272           dest -= stride[n] * extent[n];
 273           n++;
 274           if (n == dim)
 275             {
 276               dest = NULL;
 277               break;
 278             }
 279           else
 280             {
 281               count[n]++;
 282               dest += stride[n];
 283             }
 284         }
 285     }
 286 }
 287
 288 /* REAL(KIND=8) array.  */
 289
 290 #define arandom_r8 prefix(arandom_r8)
 291 void
 292 arandom_r8 (gfc_array_r8 * harv)
 293 {
 294   index_type count[GFC_MAX_DIMENSIONS - 1];
 295   index_type extent[GFC_MAX_DIMENSIONS - 1];
 296   index_type stride[GFC_MAX_DIMENSIONS - 1];
 297   index_type stride0;
 298   index_type dim;
 299   GFC_REAL_8 *dest;
 300   int n;
 301
 302   dest = harv->data;
 303
 304   if (harv->dim[0].stride == 0)
 305     harv->dim[0].stride = 1;
 306
 307   dim = GFC_DESCRIPTOR_RANK (harv);
 308
 309   for (n = 0; n < dim; n++)
 310     {
 311       count[n] = 0;
 312       stride[n] = harv->dim[n].stride;
 313       extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
 314       if (extent[n] <= 0)
 315         return;
 316     }
 317
 318   stride0 = stride[0];
 319
 320   while (dest)
 321     {
 322       /* Set the elements.  */
 323
 324       /* regenerate if we need to, may waste one 32-bit value */
 325       if ((i + 1) >= N)
 326         random_generate ();
 327
 328       /* Convert two uint32 to a REAL(KIND=8).  */
 329       *dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
 330               (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
 331       i += 2;
 332
 333       /* Advance to the next element.  */
 334       dest += stride0;
 335       count[0]++;
 336       /* Advance to the next source element.  */
 337       n = 0;
 338       while (count[n] == extent[n])
 339         {
 340           /* When we get to the end of a dimension,
 341              reset it and increment
 342              the next dimension.  */
 343           count[n] = 0;
 344           /* We could precalculate these products,
 345              but this is a less
 346              frequently used path so proabably not worth it.  */
 347           dest -= stride[n] * extent[n];
 348           n++;
 349           if (n == dim)
 350             {
 351               dest = NULL;
 352               break;
 353             }
 354           else
 355             {
 356               count[n]++;
 357               dest += stride[n];
 358             }
 359         }
 360     }
 361 }
 362