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     }
 466   while (*x == 1.0);
 467
 468 }
 469
 470 /*  This function produces a REAL(8) value from the uniform distribution
 471     with range [0,1).  */
 472
 473 void
 474 prefix(random_r8) (GFC_REAL_8 *x)
 475 {
 476
 477   GFC_UINTEGER_8 kiss;
 478
 479   do
 480     {
 481       kiss = (((GFC_UINTEGER_8)kiss_random_kernel ()) << 32)
 482              + kiss_random_kernel ();
 483       *x = (GFC_REAL_8)kiss / (GFC_REAL_8)(~(GFC_UINTEGER_8) 0);
 484     }
 485   while (*x != 0);
 486
 487 }
 488
 489 /*  This function fills a REAL(4) array with values from the uniform
 490     distribution with range [0,1).  */
 491
 492 void
 493 prefix(arandom_r4) (gfc_array_r4 *x)
 494 {
 495
 496   index_type count[GFC_MAX_DIMENSIONS - 1];
 497   index_type extent[GFC_MAX_DIMENSIONS - 1];
 498   index_type stride[GFC_MAX_DIMENSIONS - 1];
 499   index_type stride0;
 500   index_type dim;
 501   GFC_REAL_4 *dest;
 502   int n;
 503
 504   dest = x->data;
 505
 506   if (x->dim[0].stride == 0)
 507     x->dim[0].stride = 1;
 508
 509   dim = GFC_DESCRIPTOR_RANK (x);
 510
 511   for (n = 0; n < dim; n++)
 512     {
 513       count[n] = 0;
 514       stride[n] = x->dim[n].stride;
 515       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 516       if (extent[n] <= 0)
 517         return;
 518     }
 519
 520   stride0 = stride[0];
 521
 522   while (dest)
 523     {
 524       prefix(random_r4) (dest);
 525
 526       /* Advance to the next element.  */
 527       dest += stride0;
 528       count[0]++;
 529       /* Advance to the next source element.  */
 530       n = 0;
 531       while (count[n] == extent[n])
 532         {
 533           /* When we get to the end of a dimension, reset it and increment
 534              the next dimension.  */
 535           count[n] = 0;
 536           /* We could precalculate these products, but this is a less
 537              frequently used path so probably not worth it.  */
 538           dest -= stride[n] * extent[n];
 539           n++;
 540           if (n == dim)
 541             {
 542               dest = NULL;
 543               break;
 544             }
 545           else
 546             {
 547               count[n]++;
 548               dest += stride[n];
 549             }
 550         }
 551     }
 552 }
 553
 554 /*  This function fills a REAL(8) array with valuse from the uniform
 555     distribution with range [0,1).  */
 556
 557 void
 558 prefix(arandom_r8) (gfc_array_r8 *x)
 559 {
 560
 561   index_type count[GFC_MAX_DIMENSIONS - 1];
 562   index_type extent[GFC_MAX_DIMENSIONS - 1];
 563   index_type stride[GFC_MAX_DIMENSIONS - 1];
 564   index_type stride0;
 565   index_type dim;
 566   GFC_REAL_8 *dest;
 567   int n;
 568
 569   dest = x->data;
 570
 571   if (x->dim[0].stride == 0)
 572     x->dim[0].stride = 1;
 573
 574   dim = GFC_DESCRIPTOR_RANK (x);
 575
 576   for (n = 0; n < dim; n++)
 577     {
 578       count[n] = 0;
 579       stride[n] = x->dim[n].stride;
 580       extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
 581       if (extent[n] <= 0)
 582         return;
 583     }
 584
 585   stride0 = stride[0];
 586
 587   while (dest)
 588     {
 589       prefix(random_r8) (dest);
 590
 591       /* Advance to the next element.  */
 592       dest += stride0;
 593       count[0]++;
 594       /* Advance to the next source element.  */
 595       n = 0;
 596       while (count[n] == extent[n])
 597         {
 598           /* When we get to the end of a dimension, reset it and increment
 599              the next dimension.  */
 600           count[n] = 0;
 601           /* We could precalculate these products, but this is a less
 602              frequently used path so probably not worth it.  */
 603           dest -= stride[n] * extent[n];
 604           n++;
 605           if (n == dim)
 606             {
 607               dest = NULL;
 608               break;
 609             }
 610           else
 611             {
 612               count[n]++;
 613               dest += stride[n];
 614             }
 615         }
 616     }
 617 }
 618
 619 /* prefix(random_seed) is used to seed the PRNG with either a default
 620    set of seeds or user specified set of seed.  prefix(random_seed)
 621    must be called with no argument or exactly one argument.  */
 622
 623 void
 624 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 * put,
 625                      gfc_array_i4 * get)
 626 {
 627
 628   int i;
 629
 630   if (size == NULL && put == NULL && get == NULL)
 631     {
 632       /* From the standard: "If no argument is present, the processor assigns
 633          a processor-dependent value to the seed." */
 634       kiss_seed[0] = kiss_default_seed[0];
 635       kiss_seed[1] = kiss_default_seed[1];
 636       kiss_seed[2] = kiss_default_seed[2];
 637       kiss_seed[3] = kiss_default_seed[3];
 638     }
 639
 640   if (size != NULL)
 641     *size = kiss_size;
 642
 643   if (put != NULL)
 644     {
 645       /* If the rank of the array is not 1, abort */
 646       if (GFC_DESCRIPTOR_RANK (put) != 1)
 647         runtime_error ("Array rank of PUT is not 1.");
 648
 649       /* If the array is too small, abort */
 650       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
 651         runtime_error ("Array size of PUT is too small.");
 652
 653       if (put->dim[0].stride == 0)
 654         put->dim[0].stride = 1;
 655
 656       /*  This code now should do correct strides. */
 657       for (i = 0; i < kiss_size; i++)
 658         kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
 659     }
 660
 661   /* Return the seed to GET data */
 662   if (get != NULL)
 663     {
 664       /* If the rank of the array is not 1, abort. */
 665       if (GFC_DESCRIPTOR_RANK (get) != 1)
 666         runtime_error ("Array rank of GET is not 1.");
 667
 668       /* If the array is too small, abort. */
 669       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
 670         runtime_error ("Array size of GET is too small.");
 671
 672       if (get->dim[0].stride == 0)
 673         get->dim[0].stride = 1;
 674
 675       /*  This code now should do correct strides. */
 676       for (i = 0; i < kiss_size; i++)
 677         get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
 678     }
 679 }
 680
 681