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., 59 Temple Place - Suite 330,
  30 Boston, MA 02111-1307, USA.  */
  31
  32 #include "libgfortran.h"
  33
  34 extern void random_r4 (GFC_REAL_4 *);
  35 iexport_proto(random_r4);
  36
  37 extern void random_r8 (GFC_REAL_8 *);
  38 iexport_proto(random_r8);
  39
  40 extern void arandom_r4 (gfc_array_r4 *);
  41 export_proto(arandom_r4);
  42
  43 extern void arandom_r8 (gfc_array_r8 *);
  44 export_proto(arandom_r8);
  45
  46 #if 0
  47
  48 /*  The Mersenne Twister code is currently commented out due to
  49
  50     (1) Simple user specified seeds lead to really bad sequences for
  51         nearly 100000 random numbers.
  52     (2) open(), read(), and close() are not properly declared via header
  53         files.
  54     (3) The global index i is abused and causes unexpected behavior with
  55         GET and PUT.
  56     (4) See PR 15619.
  57
  58   The algorithm was taken from the paper :
  59
  60         Mersenne Twister:       623-dimensionally equidistributed
  61                                 uniform pseudorandom generator.
  62
  63         by:     Makoto Matsumoto
  64                 Takuji Nishimura
  65
  66         Which appeared in the: ACM Transactions on Modelling and Computer
  67         Simulations: Special Issue on Uniform Random Number
  68         Generation. ( Early in 1998 ).  */
  69
  70
  71 #include <stdio.h>
  72 #include <stdlib.h>
  73 #include <sys/types.h>
  74 #include <sys/stat.h>
  75 #include <fcntl.h>
  76
  77 #ifdef HAVE_UNISTD_H
  78 #include <unistd.h>
  79 #endif
  80
  81 /*Use the 'big' generator by default ( period -> 2**19937 ).  */
  82
  83 #define MT19937
  84
  85 /* Define the necessary constants for the algorithm.  */
  86
  87 #ifdef  MT19937
  88 enum constants
  89 {
  90   N = 624, M = 397, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
  91 };
  92 #define M_A     0x9908B0DF
  93 #define T_B     0x9D2C5680
  94 #define T_C     0xEFC60000
  95 #else
  96 enum constants
  97 {
  98   N = 351, M = 175, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
  99 };
 100 #define M_A     0xE4BD75F5
 101 #define T_B     0x655E5280
 102 #define T_C     0xFFD58000
 103 #endif
 104
 105 static int i = N;
 106 static unsigned int seed[N];
 107
 108 /* This is the routine which handles the seeding of the generator,
 109    and also reading and writing of the seed.  */
 110
 111 void
 112 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 113 {
 114   /* Initialize the seed in system dependent manner.  */
 115   if (get == NULL && put == NULL && size == NULL)
 116     {
 117       int fd;
 118       fd = open ("/dev/urandom", O_RDONLY);
 119       if (fd == 0)
 120         {
 121           /* We dont have urandom.  */
 122           GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed;
 123           for (i = 0; i < N; i++)
 124             {
 125               s = s * 29943829 - 1;
 126               seed[i] = s;
 127             }
 128         }
 129       else
 130         {
 131           /* Using urandom, might have a length issue.  */
 132           read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N);
 133           close (fd);
 134         }
 135       return;
 136     }
 137
 138   /* Return the size of the seed */
 139   if (size != NULL)
 140     {
 141       *size = N;
 142       return;
 143     }
 144
 145   /* if we have gotten to this pount we have a get or put
 146    * now we check it the array fulfills the demands in the standard .
 147    */
 148
 149   /* Set the seed to PUT data */
 150   if (put != NULL)
 151     {
 152       /* if the rank of the array is not 1 abort */
 153       if (GFC_DESCRIPTOR_RANK (put) != 1)
 154         abort ();
 155
 156       /* if the array is too small abort */
 157       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N)
 158         abort ();
 159
 160       /* If this is the case the array is a temporary */
 161       if (put->dim[0].stride == 0)
 162         return;
 163
 164       /*  This code now should do correct strides. */
 165       for (i = 0; i < N; i++)
 166         seed[i] = put->data[i * put->dim[0].stride];
 167     }
 168
 169   /* Return the seed to GET data */
 170   if (get != NULL)
 171     {
 172       /* if the rank of the array is not 1 abort */
 173       if (GFC_DESCRIPTOR_RANK (get) != 1)
 174         abort ();
 175
 176       /* if the array is too small abort */
 177       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N)
 178         abort ();
 179
 180       /* If this is the case the array is a temporary */
 181       if (get->dim[0].stride == 0)
 182         return;
 183
 184       /*  This code now should do correct strides. */
 185       for (i = 0; i < N; i++)
 186         get->data[i * get->dim[0].stride] = seed[i];
 187     }
 188 }
 189 iexport(random_seed);
 190
 191 /* Here is the internal routine which generates the random numbers
 192    in 'batches' based upon the need for a new batch.
 193    It's an integer based routine known as 'Mersenne Twister'.
 194    This implementation still lacks 'tempering' and a good verification,
 195    but gives very good metrics.  */
 196
 197 static void
 198 random_generate (void)
 199 {
 200   /* 32 bits.  */
 201   GFC_UINTEGER_4 y;
 202
 203   /* Generate batch of N.  */
 204   int k, m;
 205   for (k = 0, m = M; k < N - 1; k++)
 206     {
 207       y = (seed[k] & (-1 << R)) | (seed[k + 1] & ((1u << R) - 1));
 208       seed[k] = seed[m] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
 209       if (++m >= N)
 210         m = 0;
 211     }
 212
 213   y = (seed[N - 1] & (-1 << R)) | (seed[0] & ((1u << R) - 1));
 214   seed[N - 1] = seed[M - 1] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
 215   i = 0;
 216 }
 217
 218 /* A routine to return a REAL(KIND=4).  */
 219
 220 void
 221 random_r4 (GFC_REAL_4 * harv)
 222 {
 223   /* Regenerate if we need to.  */
 224   if (i >= N)
 225     random_generate ();
 226
 227   /* Convert uint32 to REAL(KIND=4).  */
 228   *harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
 229                         (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
 230 }
 231 iexport(random_r4);
 232
 233 /* A routine to return a REAL(KIND=8).  */
 234
 235 void
 236 random_r8 (GFC_REAL_8 * harv)
 237 {
 238   /* Regenerate if we need to, may waste one 32-bit value.  */
 239   if ((i + 1) >= N)
 240     random_generate ();
 241
 242   /* Convert two uint32 to a REAL(KIND=8).  */
 243   *harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
 244           (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
 245   i += 2;
 246 }
 247 iexport(random_r8);
 248
 249 /* Code to handle arrays will follow here.  */
 250
 251 /* REAL(KIND=4) REAL array.  */
 252
 253 void
 254 arandom_r4 (gfc_array_r4 * harv)
 255 {
 256   index_type count[GFC_MAX_DIMENSIONS];
 257   index_type extent[GFC_MAX_DIMENSIONS];
 258   index_type stride[GFC_MAX_DIMENSIONS];
 259   index_type stride0;
 260   index_type dim;
 261   GFC_REAL_4 *dest;
 262   int n;
 263
 264   dest = harv->data;
 265
 266   if (harv->dim[0].stride == 0)
 267     harv->dim[0].stride = 1;
 268
 269   dim = GFC_DESCRIPTOR_RANK (harv);
 270
 271   for (n = 0; n < dim; n++)
 272     {
 273       count[n] = 0;
 274       stride[n] = harv->dim[n].stride;
 275       extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
 276       if (extent[n] <= 0)
 277         return;
 278     }
 279
 280   stride0 = stride[0];
 281
 282   while (dest)
 283     {
 284       /* Set the elements.  */
 285
 286       /* regenerate if we need to */
 287       if (i >= N)
 288         random_generate ();
 289
 290       /* Convert uint32 to float in a hopefully g95 compiant manner */
 291       *dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
 292                             (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
 293
 294       /* Advance to the next element.  */
 295       dest += stride0;
 296       count[0]++;
 297       /* Advance to the next source element.  */
 298       n = 0;
 299       while (count[n] == extent[n])
 300         {
 301           /* When we get to the end of a dimension,
 302              reset it and increment
 303              the next dimension.  */
 304           count[n] = 0;
 305           /* We could precalculate these products,
 306              but this is a less
 307              frequently used path so proabably not worth it.  */
 308           dest -= stride[n] * extent[n];
 309           n++;
 310           if (n == dim)
 311             {
 312               dest = NULL;
 313               break;
 314             }
 315           else
 316             {
 317               count[n]++;
 318               dest += stride[n];
 319             }
 320         }
 321     }
 322 }
 323
 324 /* REAL(KIND=8) array.  */
 325
 326 void
 327 arandom_r8 (gfc_array_r8 * harv)
 328 {
 329   index_type count[GFC_MAX_DIMENSIONS];
 330   index_type extent[GFC_MAX_DIMENSIONS];
 331   index_type stride[GFC_MAX_DIMENSIONS];
 332   index_type stride0;
 333   index_type dim;
 334   GFC_REAL_8 *dest;
 335   int n;
 336
 337   dest = harv->data;
 338
 339   if (harv->dim[0].stride == 0)
 340     harv->dim[0].stride = 1;
 341
 342   dim = GFC_DESCRIPTOR_RANK (harv);
 343
 344   for (n = 0; n < dim; n++)
 345     {
 346       count[n] = 0;
 347       stride[n] = harv->dim[n].stride;
 348       extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
 349       if (extent[n] <= 0)
 350         return;
 351     }
 352
 353   stride0 = stride[0];
 354
 355   while (dest)
 356     {
 357       /* Set the elements.  */
 358
 359       /* regenerate if we need to, may waste one 32-bit value */
 360       if ((i + 1) >= N)
 361         random_generate ();
 362
 363       /* Convert two uint32 to a REAL(KIND=8).  */
 364       *dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
 365               (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
 366       i += 2;
 367
 368       /* Advance to the next element.  */
 369       dest += stride0;
 370       count[0]++;
 371       /* Advance to the next source element.  */
 372       n = 0;
 373       while (count[n] == extent[n])
 374         {
 375           /* When we get to the end of a dimension,
 376              reset it and increment
 377              the next dimension.  */
 378           count[n] = 0;
 379           /* We could precalculate these products,
 380              but this is a less
 381              frequently used path so proabably not worth it.  */
 382           dest -= stride[n] * extent[n];
 383           n++;
 384           if (n == dim)
 385             {
 386               dest = NULL;
 387               break;
 388             }
 389           else
 390             {
 391               count[n]++;
 392               dest += stride[n];
 393             }
 394         }
 395     }
 396 }
 397
 398 #else
 399
 400 /* George Marsaglia's KISS (Keep It Simple Stupid) random number generator.
 401
 402    This PRNG combines:
 403
 404    (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
 405        of 2^32,
 406    (2) A 3-shift shift-register generator with a period of 2^32-1,
 407    (3) Two 16-bit multiply-with-carry generators with a period of
 408        597273182964842497 > 2^59.
 409
 410    The overall period exceeds 2^123.
 411
 412    http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
 413
 414    The above web site has an archive of a newsgroup posting from George
 415    Marsaglia with the statement:
 416
 417    Subject: Random numbers for C: Improvements.
 418    Date: Fri, 15 Jan 1999 11:41:47 -0500
 419    From: George Marsaglia <geo@stat.fsu.edu>
 420    Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
 421    References: <369B5E30.65A55FD1@stat.fsu.edu>
 422    Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
 423    Lines: 93
 424
 425    As I hoped, several suggestions have led to
 426    improvements in the code for RNG's I proposed for
 427    use in C. (See the thread "Random numbers for C: Some
 428    suggestions" in previous postings.) The improved code
 429    is listed below.
 430
 431    A question of copyright has also been raised.  Unlike
 432    DIEHARD, there is no copyright on the code below. You
 433    are free to use it in any way you want, but you may
 434    wish to acknowledge the source, as a courtesy.
 435
 436 "There is no copyright on the code below." included the original
 437 KISS algorithm.  */
 438
 439 #define GFC_SL(k, n)    ((k)^((k)<<(n)))
 440 #define GFC_SR(k, n)    ((k)^((k)>>(n)))
 441
 442 static const GFC_INTEGER_4 kiss_size = 4;
 443 #define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069}
 444 static const GFC_UINTEGER_4 kiss_default_seed[4] = KISS_DEFAULT_SEED;
 445 static GFC_UINTEGER_4 kiss_seed[4] = KISS_DEFAULT_SEED;
 446
 447 /* kiss_random_kernel() returns an integer value in the range of
 448    (0, GFC_UINTEGER_4_HUGE].  The distribution of pseudorandom numbers
 449    should be uniform.  */
 450
 451 static GFC_UINTEGER_4
 452 kiss_random_kernel(void)
 453 {
 454   GFC_UINTEGER_4 kiss;
 455
 456   kiss_seed[0] = 69069 * kiss_seed[0] + 1327217885;
 457   kiss_seed[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed[1],13),17),5);
 458   kiss_seed[2] = 18000 * (kiss_seed[2] & 65535) + (kiss_seed[2] >> 16);
 459   kiss_seed[3] = 30903 * (kiss_seed[3] & 65535) + (kiss_seed[3] >> 16);
 460   kiss = kiss_seed[0] + kiss_seed[1] + (kiss_seed[2] << 16) + kiss_seed[3];
 461
 462   return kiss;
 463 }
 464
 465 /*  This function produces a REAL(4) value from the uniform distribution
 466     with range [0,1).  */
 467
 468 void
 469 random_r4 (GFC_REAL_4 *x)
 470 {
 471   GFC_UINTEGER_4 kiss;
 472
 473   kiss = kiss_random_kernel ();
 474   /* Burn a random number, so the REAL*4 and REAL*8 functions
 475      produce similar sequences of random numbers.  */
 476   kiss_random_kernel ();
 477   *x = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
 478 }
 479 iexport(random_r4);
 480
 481 /*  This function produces a REAL(8) value from the uniform distribution
 482     with range [0,1).  */
 483
 484 void
 485 random_r8 (GFC_REAL_8 *x)
 486 {
 487   GFC_UINTEGER_8 kiss;
 488
 489   kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
 490   kiss += kiss_random_kernel ();
 491   *x = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
 492 }
 493 iexport(random_r8);
 494
 495 /*  This function fills a REAL(4) array with values from the uniform
 496     distribution with range [0,1).  */
 497
 498 void
 499 arandom_r4 (gfc_array_r4 *x)
 500 {
 501   index_type count[GFC_MAX_DIMENSIONS];
 502   index_type extent[GFC_MAX_DIMENSIONS];
 503   index_type stride[GFC_MAX_DIMENSIONS];
 504   index_type stride0;
 505   index_type dim;
 506   GFC_REAL_4 *dest;
 507   int n;
 508
 509   dest = x->data;
 510
 511   if (x->dim[0].stride == 0)
 512     x->dim[0].stride = 1;
 513
 514   dim = GFC_DESCRIPTOR_RANK (x);
 515
 516   for (n = 0; n < dim; n++)
 517     {
 518       count[n] = 0;
 519       stride[n] = x->dim[n].stride;
 520       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 521       if (extent[n] <= 0)
 522         return;
 523     }
 524
 525   stride0 = stride[0];
 526
 527   while (dest)
 528     {
 529       random_r4 (dest);
 530
 531       /* Advance to the next element.  */
 532       dest += stride0;
 533       count[0]++;
 534       /* Advance to the next source element.  */
 535       n = 0;
 536       while (count[n] == extent[n])
 537         {
 538           /* When we get to the end of a dimension, reset it and increment
 539              the next dimension.  */
 540           count[n] = 0;
 541           /* We could precalculate these products, but this is a less
 542              frequently used path so probably not worth it.  */
 543           dest -= stride[n] * extent[n];
 544           n++;
 545           if (n == dim)
 546             {
 547               dest = NULL;
 548               break;
 549             }
 550           else
 551             {
 552               count[n]++;
 553               dest += stride[n];
 554             }
 555         }
 556     }
 557 }
 558
 559 /*  This function fills a REAL(8) array with values from the uniform
 560     distribution with range [0,1).  */
 561
 562 void
 563 arandom_r8 (gfc_array_r8 *x)
 564 {
 565   index_type count[GFC_MAX_DIMENSIONS];
 566   index_type extent[GFC_MAX_DIMENSIONS];
 567   index_type stride[GFC_MAX_DIMENSIONS];
 568   index_type stride0;
 569   index_type dim;
 570   GFC_REAL_8 *dest;
 571   int n;
 572
 573   dest = x->data;
 574
 575   if (x->dim[0].stride == 0)
 576     x->dim[0].stride = 1;
 577
 578   dim = GFC_DESCRIPTOR_RANK (x);
 579
 580   for (n = 0; n < dim; n++)
 581     {
 582       count[n] = 0;
 583       stride[n] = x->dim[n].stride;
 584       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 585       if (extent[n] <= 0)
 586         return;
 587     }
 588
 589   stride0 = stride[0];
 590
 591   while (dest)
 592     {
 593       random_r8 (dest);
 594
 595       /* Advance to the next element.  */
 596       dest += stride0;
 597       count[0]++;
 598       /* Advance to the next source element.  */
 599       n = 0;
 600       while (count[n] == extent[n])
 601         {
 602           /* When we get to the end of a dimension, reset it and increment
 603              the next dimension.  */
 604           count[n] = 0;
 605           /* We could precalculate these products, but this is a less
 606              frequently used path so probably not worth it.  */
 607           dest -= stride[n] * extent[n];
 608           n++;
 609           if (n == dim)
 610             {
 611               dest = NULL;
 612               break;
 613             }
 614           else
 615             {
 616               count[n]++;
 617               dest += stride[n];
 618             }
 619         }
 620     }
 621 }
 622
 623 /* random_seed is used to seed the PRNG with either a default
 624    set of seeds or user specified set of seeds.  random_seed
 625    must be called with no argument or exactly one argument.  */
 626
 627 void
 628 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 629 {
 630   int i;
 631
 632   if (size == NULL && put == NULL && get == NULL)
 633     {
 634       /* From the standard: "If no argument is present, the processor assigns
 635          a processor-dependent value to the seed."  */
 636       kiss_seed[0] = kiss_default_seed[0];
 637       kiss_seed[1] = kiss_default_seed[1];
 638       kiss_seed[2] = kiss_default_seed[2];
 639       kiss_seed[3] = kiss_default_seed[3];
 640     }
 641
 642   if (size != NULL)
 643     *size = kiss_size;
 644
 645   if (put != NULL)
 646     {
 647       /* If the rank of the array is not 1, abort.  */
 648       if (GFC_DESCRIPTOR_RANK (put) != 1)
 649         runtime_error ("Array rank of PUT is not 1.");
 650
 651       /* If the array is too small, abort.  */
 652       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
 653         runtime_error ("Array size of PUT is too small.");
 654
 655       if (put->dim[0].stride == 0)
 656         put->dim[0].stride = 1;
 657
 658       /*  This code now should do correct strides.  */
 659       for (i = 0; i < kiss_size; i++)
 660         kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
 661     }
 662
 663   /* Return the seed to GET data.  */
 664   if (get != NULL)
 665     {
 666       /* If the rank of the array is not 1, abort.  */
 667       if (GFC_DESCRIPTOR_RANK (get) != 1)
 668         runtime_error ("Array rank of GET is not 1.");
 669
 670       /* If the array is too small, abort.  */
 671       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
 672         runtime_error ("Array size of GET is too small.");
 673
 674       if (get->dim[0].stride == 0)
 675         get->dim[0].stride = 1;
 676
 677       /*  This code now should do correct strides.  */
 678       for (i = 0; i < kiss_size; i++)
 679         get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
 680     }
 681 }
 682 iexport(random_seed);
 683
 684 #endif /* mersenne twister */