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