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