libgfortran/intrinsics/random.c

   1 /* Implementation of the RANDOM intrinsics
   2    Copyright 2002, 2004, 2005, 2006 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   dim = GFC_DESCRIPTOR_RANK (x);
 192
 193   for (n = 0; n < dim; n++)
 194     {
 195       count[n] = 0;
 196       stride[n] = x->dim[n].stride;
 197       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 198       if (extent[n] <= 0)
 199         return;
 200     }
 201
 202   stride0 = stride[0];
 203
 204   __gthread_mutex_lock (&random_lock);
 205
 206   while (dest)
 207     {
 208       /* random_r4 (dest); */
 209       kiss = kiss_random_kernel ();
 210       /* Burn a random number, so the REAL*4 and REAL*8 functions
 211          produce similar sequences of random numbers.  */
 212       kiss_random_kernel ();
 213       *dest = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
 214
 215       /* Advance to the next element.  */
 216       dest += stride0;
 217       count[0]++;
 218       /* Advance to the next source element.  */
 219       n = 0;
 220       while (count[n] == extent[n])
 221         {
 222           /* When we get to the end of a dimension, reset it and increment
 223              the next dimension.  */
 224           count[n] = 0;
 225           /* We could precalculate these products, but this is a less
 226              frequently used path so probably not worth it.  */
 227           dest -= stride[n] * extent[n];
 228           n++;
 229           if (n == dim)
 230             {
 231               dest = NULL;
 232               break;
 233             }
 234           else
 235             {
 236               count[n]++;
 237               dest += stride[n];
 238             }
 239         }
 240     }
 241   __gthread_mutex_unlock (&random_lock);
 242 }
 243
 244 /*  This function fills a REAL(8) array with values from the uniform
 245     distribution with range [0,1).  */
 246
 247 void
 248 arandom_r8 (gfc_array_r8 *x)
 249 {
 250   index_type count[GFC_MAX_DIMENSIONS];
 251   index_type extent[GFC_MAX_DIMENSIONS];
 252   index_type stride[GFC_MAX_DIMENSIONS];
 253   index_type stride0;
 254   index_type dim;
 255   GFC_REAL_8 *dest;
 256   GFC_UINTEGER_8 kiss;
 257   int n;
 258
 259   dest = x->data;
 260
 261   dim = GFC_DESCRIPTOR_RANK (x);
 262
 263   for (n = 0; n < dim; n++)
 264     {
 265       count[n] = 0;
 266       stride[n] = x->dim[n].stride;
 267       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 268       if (extent[n] <= 0)
 269         return;
 270     }
 271
 272   stride0 = stride[0];
 273
 274   __gthread_mutex_lock (&random_lock);
 275
 276   while (dest)
 277     {
 278       /* random_r8 (dest); */
 279       kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
 280       kiss += kiss_random_kernel ();
 281       *dest = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
 282
 283       /* Advance to the next element.  */
 284       dest += stride0;
 285       count[0]++;
 286       /* Advance to the next source element.  */
 287       n = 0;
 288       while (count[n] == extent[n])
 289         {
 290           /* When we get to the end of a dimension, reset it and increment
 291              the next dimension.  */
 292           count[n] = 0;
 293           /* We could precalculate these products, but this is a less
 294              frequently used path so probably not worth it.  */
 295           dest -= stride[n] * extent[n];
 296           n++;
 297           if (n == dim)
 298             {
 299               dest = NULL;
 300               break;
 301             }
 302           else
 303             {
 304               count[n]++;
 305               dest += stride[n];
 306             }
 307         }
 308     }
 309   __gthread_mutex_unlock (&random_lock);
 310 }
 311
 312 /* random_seed is used to seed the PRNG with either a default
 313    set of seeds or user specified set of seeds.  random_seed
 314    must be called with no argument or exactly one argument.  */
 315
 316 void
 317 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 318 {
 319   int i;
 320
 321   __gthread_mutex_lock (&random_lock);
 322
 323   if (size == NULL && put == NULL && get == NULL)
 324     {
 325       /* From the standard: "If no argument is present, the processor assigns
 326          a processor-dependent value to the seed."  */
 327       kiss_seed[0] = kiss_default_seed[0];
 328       kiss_seed[1] = kiss_default_seed[1];
 329       kiss_seed[2] = kiss_default_seed[2];
 330       kiss_seed[3] = kiss_default_seed[3];
 331     }
 332
 333   if (size != NULL)
 334     *size = kiss_size;
 335
 336   if (put != NULL)
 337     {
 338       /* If the rank of the array is not 1, abort.  */
 339       if (GFC_DESCRIPTOR_RANK (put) != 1)
 340         runtime_error ("Array rank of PUT is not 1.");
 341
 342       /* If the array is too small, abort.  */
 343       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
 344         runtime_error ("Array size of PUT is too small.");
 345
 346       /*  This code now should do correct strides.  */
 347       for (i = 0; i < kiss_size; i++)
 348         kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
 349     }
 350
 351   /* Return the seed to GET data.  */
 352   if (get != NULL)
 353     {
 354       /* If the rank of the array is not 1, abort.  */
 355       if (GFC_DESCRIPTOR_RANK (get) != 1)
 356         runtime_error ("Array rank of GET is not 1.");
 357
 358       /* If the array is too small, abort.  */
 359       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
 360         runtime_error ("Array size of GET is too small.");
 361
 362       /*  This code now should do correct strides.  */
 363       for (i = 0; i < kiss_size; i++)
 364         get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
 365     }
 366
 367   __gthread_mutex_unlock (&random_lock);
 368 }
 369 iexport(random_seed);
 370
 371
 372 #ifndef __GTHREAD_MUTEX_INIT
 373 static void __attribute__((constructor))
 374 init (void)
 375 {
 376   __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
 377 }
 378 #endif