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