old-autovect-branch/libgfortran/intrinsics/random.c

   1 /* Implementation of the RANDOM intrinsics
   2    Copyright 2002, 2004, 2005 Free Software Foundation, Inc.
   3    Contributed by Lars Segerlund <seger@linuxmail.org>
   4    and Steve Kargl.
   5
   6 This file is part of the GNU Fortran 95 runtime library (libgfortran).
   7
   8 Libgfortran is free software; you can redistribute it and/or
   9 modify it under the terms of the GNU General Public
  10 License as published by the Free Software Foundation; either
  11 version 2 of the License, or (at your option) any later version.
  12
  13 In addition to the permissions in the GNU General Public License, the
  14 Free Software Foundation gives you unlimited permission to link the
  15 compiled version of this file into combinations with other programs,
  16 and to distribute those combinations without any restriction coming
  17 from the use of this file.  (The General Public License restrictions
  18 do apply in other respects; for example, they cover modification of
  19 the file, and distribution when not linked into a combine
  20 executable.)
  21
  22 Ligbfortran is distributed in the hope that it will be useful,
  23 but WITHOUT ANY WARRANTY; without even the implied warranty of
  24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  25 GNU General Public License for more details.
  26
  27 You should have received a copy of the GNU General Public
  28 License along with libgfortran; see the file COPYING.  If not,
  29 write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  30 Boston, MA 02110-1301, USA.  */
  31
  32 #include "config.h"
  33 #include "libgfortran.h"
  34 #include "../io/io.h"
  35
  36 extern void random_r4 (GFC_REAL_4 *);
  37 iexport_proto(random_r4);
  38
  39 extern void random_r8 (GFC_REAL_8 *);
  40 iexport_proto(random_r8);
  41
  42 extern void arandom_r4 (gfc_array_r4 *);
  43 export_proto(arandom_r4);
  44
  45 extern void arandom_r8 (gfc_array_r8 *);
  46 export_proto(arandom_r8);
  47
  48 #ifdef __GTHREAD_MUTEX_INIT
  49 static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
  50 #else
  51 static __gthread_mutex_t random_lock;
  52 #endif
  53
  54
  55 /* libgfortran previously had a Mersenne Twister, taken from the paper:
  56
  57         Mersenne Twister:       623-dimensionally equidistributed
  58                                 uniform pseudorandom generator.
  59
  60         by Makoto Matsumoto & Takuji Nishimura
  61         which appeared in the: ACM Transactions on Modelling and Computer
  62         Simulations: Special Issue on Uniform Random Number
  63         Generation. ( Early in 1998 ).
  64
  65    The Mersenne Twister code was replaced due to
  66
  67     (1) Simple user specified seeds lead to really bad sequences for
  68         nearly 100000 random numbers.
  69     (2) open(), read(), and close() were not properly declared via header
  70         files.
  71     (3) The global index i was abused and caused unexpected behavior with
  72         GET and PUT.
  73     (4) See PR 15619.
  74
  75
  76    libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
  77    random number generator.  This PRNG combines:
  78
  79    (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
  80        of 2^32,
  81    (2) A 3-shift shift-register generator with a period of 2^32-1,
  82    (3) Two 16-bit multiply-with-carry generators with a period of
  83        597273182964842497 > 2^59.
  84
  85    The overall period exceeds 2^123.
  86
  87    http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
  88
  89    The above web site has an archive of a newsgroup posting from George
  90    Marsaglia with the statement:
  91
  92    Subject: Random numbers for C: Improvements.
  93    Date: Fri, 15 Jan 1999 11:41:47 -0500
  94    From: George Marsaglia <geo@stat.fsu.edu>
  95    Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
  96    References: <369B5E30.65A55FD1@stat.fsu.edu>
  97    Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
  98    Lines: 93
  99
 100    As I hoped, several suggestions have led to
 101    improvements in the code for RNG's I proposed for
 102    use in C. (See the thread "Random numbers for C: Some
 103    suggestions" in previous postings.) The improved code
 104    is listed below.
 105
 106    A question of copyright has also been raised.  Unlike
 107    DIEHARD, there is no copyright on the code below. You
 108    are free to use it in any way you want, but you may
 109    wish to acknowledge the source, as a courtesy.
 110
 111 "There is no copyright on the code below." included the original
 112 KISS algorithm.  */
 113
 114 #define GFC_SL(k, n)    ((k)^((k)<<(n)))
 115 #define GFC_SR(k, n)    ((k)^((k)>>(n)))
 116
 117 static const GFC_INTEGER_4 kiss_size = 4;
 118 #define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069}
 119 static const GFC_UINTEGER_4 kiss_default_seed[4] = KISS_DEFAULT_SEED;
 120 static GFC_UINTEGER_4 kiss_seed[4] = KISS_DEFAULT_SEED;
 121
 122 /* kiss_random_kernel() returns an integer value in the range of
 123    (0, GFC_UINTEGER_4_HUGE].  The distribution of pseudorandom numbers
 124    should be uniform.  */
 125
 126 static GFC_UINTEGER_4
 127 kiss_random_kernel(void)
 128 {
 129   GFC_UINTEGER_4 kiss;
 130
 131   kiss_seed[0] = 69069 * kiss_seed[0] + 1327217885;
 132   kiss_seed[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed[1],13),17),5);
 133   kiss_seed[2] = 18000 * (kiss_seed[2] & 65535) + (kiss_seed[2] >> 16);
 134   kiss_seed[3] = 30903 * (kiss_seed[3] & 65535) + (kiss_seed[3] >> 16);
 135   kiss = kiss_seed[0] + kiss_seed[1] + (kiss_seed[2] << 16) + kiss_seed[3];
 136
 137   return kiss;
 138 }
 139
 140 /*  This function produces a REAL(4) value from the uniform distribution
 141     with range [0,1).  */
 142
 143 void
 144 random_r4 (GFC_REAL_4 *x)
 145 {
 146   GFC_UINTEGER_4 kiss;
 147
 148   __gthread_mutex_lock (&random_lock);
 149   kiss = kiss_random_kernel ();
 150   /* Burn a random number, so the REAL*4 and REAL*8 functions
 151      produce similar sequences of random numbers.  */
 152   kiss_random_kernel ();
 153   *x = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
 154   __gthread_mutex_unlock (&random_lock);
 155 }
 156 iexport(random_r4);
 157
 158 /*  This function produces a REAL(8) value from the uniform distribution
 159     with range [0,1).  */
 160
 161 void
 162 random_r8 (GFC_REAL_8 *x)
 163 {
 164   GFC_UINTEGER_8 kiss;
 165
 166   __gthread_mutex_lock (&random_lock);
 167   kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
 168   kiss += kiss_random_kernel ();
 169   *x = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
 170   __gthread_mutex_unlock (&random_lock);
 171 }
 172 iexport(random_r8);
 173
 174 /*  This function fills a REAL(4) array with values from the uniform
 175     distribution with range [0,1).  */
 176
 177 void
 178 arandom_r4 (gfc_array_r4 *x)
 179 {
 180   index_type count[GFC_MAX_DIMENSIONS];
 181   index_type extent[GFC_MAX_DIMENSIONS];
 182   index_type stride[GFC_MAX_DIMENSIONS];
 183   index_type stride0;
 184   index_type dim;
 185   GFC_REAL_4 *dest;
 186   GFC_UINTEGER_4 kiss;
 187   int n;
 188
 189   dest = x->data;
 190
 191   if (x->dim[0].stride == 0)
 192     x->dim[0].stride = 1;
 193
 194   dim = GFC_DESCRIPTOR_RANK (x);
 195
 196   for (n = 0; n < dim; n++)
 197     {
 198       count[n] = 0;
 199       stride[n] = x->dim[n].stride;
 200       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 201       if (extent[n] <= 0)
 202         return;
 203     }
 204
 205   stride0 = stride[0];
 206
 207   __gthread_mutex_lock (&random_lock);
 208
 209   while (dest)
 210     {
 211       /* random_r4 (dest); */
 212       kiss = kiss_random_kernel ();
 213       /* Burn a random number, so the REAL*4 and REAL*8 functions
 214          produce similar sequences of random numbers.  */
 215       kiss_random_kernel ();
 216       *dest = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
 217
 218       /* Advance to the next element.  */
 219       dest += stride0;
 220       count[0]++;
 221       /* Advance to the next source element.  */
 222       n = 0;
 223       while (count[n] == extent[n])
 224         {
 225           /* When we get to the end of a dimension, reset it and increment
 226              the next dimension.  */
 227           count[n] = 0;
 228           /* We could precalculate these products, but this is a less
 229              frequently used path so probably not worth it.  */
 230           dest -= stride[n] * extent[n];
 231           n++;
 232           if (n == dim)
 233             {
 234               dest = NULL;
 235               break;
 236             }
 237           else
 238             {
 239               count[n]++;
 240               dest += stride[n];
 241             }
 242         }
 243     }
 244   __gthread_mutex_unlock (&random_lock);
 245 }
 246
 247 /*  This function fills a REAL(8) array with values from the uniform
 248     distribution with range [0,1).  */
 249
 250 void
 251 arandom_r8 (gfc_array_r8 *x)
 252 {
 253   index_type count[GFC_MAX_DIMENSIONS];
 254   index_type extent[GFC_MAX_DIMENSIONS];
 255   index_type stride[GFC_MAX_DIMENSIONS];
 256   index_type stride0;
 257   index_type dim;
 258   GFC_REAL_8 *dest;
 259   GFC_UINTEGER_8 kiss;
 260   int n;
 261
 262   dest = x->data;
 263
 264   if (x->dim[0].stride == 0)
 265     x->dim[0].stride = 1;
 266
 267   dim = GFC_DESCRIPTOR_RANK (x);
 268
 269   for (n = 0; n < dim; n++)
 270     {
 271       count[n] = 0;
 272       stride[n] = x->dim[n].stride;
 273       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 274       if (extent[n] <= 0)
 275         return;
 276     }
 277
 278   stride0 = stride[0];
 279
 280   __gthread_mutex_lock (&random_lock);
 281
 282   while (dest)
 283     {
 284       /* random_r8 (dest); */
 285       kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
 286       kiss += kiss_random_kernel ();
 287       *dest = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
 288
 289       /* Advance to the next element.  */
 290       dest += stride0;
 291       count[0]++;
 292       /* Advance to the next source element.  */
 293       n = 0;
 294       while (count[n] == extent[n])
 295         {
 296           /* When we get to the end of a dimension, reset it and increment
 297              the next dimension.  */
 298           count[n] = 0;
 299           /* We could precalculate these products, but this is a less
 300              frequently used path so probably not worth it.  */
 301           dest -= stride[n] * extent[n];
 302           n++;
 303           if (n == dim)
 304             {
 305               dest = NULL;
 306               break;
 307             }
 308           else
 309             {
 310               count[n]++;
 311               dest += stride[n];
 312             }
 313         }
 314     }
 315   __gthread_mutex_unlock (&random_lock);
 316 }
 317
 318 /* random_seed is used to seed the PRNG with either a default
 319    set of seeds or user specified set of seeds.  random_seed
 320    must be called with no argument or exactly one argument.  */
 321
 322 void
 323 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 324 {
 325   int i;
 326
 327   __gthread_mutex_lock (&random_lock);
 328
 329   if (size == NULL && put == NULL && get == NULL)
 330     {
 331       /* From the standard: "If no argument is present, the processor assigns
 332          a processor-dependent value to the seed."  */
 333       kiss_seed[0] = kiss_default_seed[0];
 334       kiss_seed[1] = kiss_default_seed[1];
 335       kiss_seed[2] = kiss_default_seed[2];
 336       kiss_seed[3] = kiss_default_seed[3];
 337     }
 338
 339   if (size != NULL)
 340     *size = kiss_size;
 341
 342   if (put != NULL)
 343     {
 344       /* If the rank of the array is not 1, abort.  */
 345       if (GFC_DESCRIPTOR_RANK (put) != 1)
 346         runtime_error ("Array rank of PUT is not 1.");
 347
 348       /* If the array is too small, abort.  */
 349       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
 350         runtime_error ("Array size of PUT is too small.");
 351
 352       if (put->dim[0].stride == 0)
 353         put->dim[0].stride = 1;
 354
 355       /*  This code now should do correct strides.  */
 356       for (i = 0; i < kiss_size; i++)
 357         kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
 358     }
 359
 360   /* Return the seed to GET data.  */
 361   if (get != NULL)
 362     {
 363       /* If the rank of the array is not 1, abort.  */
 364       if (GFC_DESCRIPTOR_RANK (get) != 1)
 365         runtime_error ("Array rank of GET is not 1.");
 366
 367       /* If the array is too small, abort.  */
 368       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
 369         runtime_error ("Array size of GET is too small.");
 370
 371       if (get->dim[0].stride == 0)
 372         get->dim[0].stride = 1;
 373
 374       /*  This code now should do correct strides.  */
 375       for (i = 0; i < kiss_size; i++)
 376         get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
 377     }
 378
 379   __gthread_mutex_unlock (&random_lock);
 380 }
 381 iexport(random_seed);
 382
 383
 384 #ifndef __GTHREAD_MUTEX_INIT
 385 static void __attribute__((constructor))
 386 init (void)
 387 {
 388   __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
 389 }
 390 #endif