libgfortran/intrinsics/random.c

   1 /* Implementation of the RANDOM intrinsics
   2    Copyright 2002, 2004, 2005, 2006, 2007, 2009, 2010
   3    Free Software Foundation, Inc.
   4    Contributed by Lars Segerlund <seger@linuxmail.org>
   5    and Steve Kargl.
   6
   7 This file is part of the GNU Fortran 95 runtime library (libgfortran).
   8
   9 Libgfortran is free software; you can redistribute it and/or
  10 modify it under the terms of the GNU General Public
  11 License as published by the Free Software Foundation; either
  12 version 3 of the License, or (at your option) any later version.
  13
  14 Ligbfortran is distributed in the hope that it will be useful,
  15 but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 GNU General Public License for more details.
  18
  19 Under Section 7 of GPL version 3, you are granted additional
  20 permissions described in the GCC Runtime Library Exception, version
  21 3.1, as published by the Free Software Foundation.
  22
  23 You should have received a copy of the GNU General Public License and
  24 a copy of the GCC Runtime Library Exception along with this program;
  25 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  26 <http://www.gnu.org/licenses/>.  */
  27
  28 #include "libgfortran.h"
  29 #include <gthr.h>
  30 #include <string.h>
  31
  32 extern void random_r4 (GFC_REAL_4 *);
  33 iexport_proto(random_r4);
  34
  35 extern void random_r8 (GFC_REAL_8 *);
  36 iexport_proto(random_r8);
  37
  38 extern void arandom_r4 (gfc_array_r4 *);
  39 export_proto(arandom_r4);
  40
  41 extern void arandom_r8 (gfc_array_r8 *);
  42 export_proto(arandom_r8);
  43
  44 #ifdef HAVE_GFC_REAL_10
  45
  46 extern void random_r10 (GFC_REAL_10 *);
  47 iexport_proto(random_r10);
  48
  49 extern void arandom_r10 (gfc_array_r10 *);
  50 export_proto(arandom_r10);
  51
  52 #endif
  53
  54 #ifdef HAVE_GFC_REAL_16
  55
  56 extern void random_r16 (GFC_REAL_16 *);
  57 iexport_proto(random_r16);
  58
  59 extern void arandom_r16 (gfc_array_r16 *);
  60 export_proto(arandom_r16);
  61
  62 #endif
  63
  64 #ifdef __GTHREAD_MUTEX_INIT
  65 static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
  66 #else
  67 static __gthread_mutex_t random_lock;
  68 #endif
  69
  70 /* Helper routines to map a GFC_UINTEGER_* to the corresponding
  71    GFC_REAL_* types in the range of [0,1).  If GFC_REAL_*_RADIX are 2
  72    or 16, respectively, we mask off the bits that don't fit into the
  73    correct GFC_REAL_*, convert to the real type, then multiply by the
  74    correct offset.  */
  75
  76
  77 static inline void
  78 rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v)
  79 {
  80   GFC_UINTEGER_4 mask;
  81 #if GFC_REAL_4_RADIX == 2
  82   mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS);
  83 #elif GFC_REAL_4_RADIX == 16
  84   mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4);
  85 #else
  86 #error "GFC_REAL_4_RADIX has unknown value"
  87 #endif
  88   v = v & mask;
  89   *f = (GFC_REAL_4) v * GFC_REAL_4_LITERAL(0x1.p-32);
  90 }
  91
  92 static inline void
  93 rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v)
  94 {
  95   GFC_UINTEGER_8 mask;
  96 #if GFC_REAL_8_RADIX == 2
  97   mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_8_DIGITS);
  98 #elif GFC_REAL_8_RADIX == 16
  99   mask = ~ (GFC_UINTEGER_8) 0u << (16 - GFC_REAL_8_DIGITS) * 4);
 100 #else
 101 #error "GFC_REAL_8_RADIX has unknown value"
 102 #endif
 103   v = v & mask;
 104   *f = (GFC_REAL_8) v * GFC_REAL_8_LITERAL(0x1.p-64);
 105 }
 106
 107 #ifdef HAVE_GFC_REAL_10
 108
 109 static inline void
 110 rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v)
 111 {
 112   GFC_UINTEGER_8 mask;
 113 #if GFC_REAL_10_RADIX == 2
 114   mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_10_DIGITS);
 115 #elif GFC_REAL_10_RADIX == 16
 116   mask = ~ (GFC_UINTEGER_10) 0u << ((16 - GFC_REAL_10_DIGITS) * 4);
 117 #else
 118 #error "GFC_REAL_10_RADIX has unknown value"
 119 #endif
 120   v = v & mask;
 121   *f = (GFC_REAL_10) v * GFC_REAL_10_LITERAL(0x1.p-64);
 122 }
 123 #endif
 124
 125 #ifdef HAVE_GFC_REAL_16
 126
 127 /* For REAL(KIND=16), we only need to mask off the lower bits.  */
 128
 129 static inline void
 130 rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
 131 {
 132   GFC_UINTEGER_8 mask;
 133 #if GFC_REAL_16_RADIX == 2
 134   mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_16_DIGITS);
 135 #elif GFC_REAL_16_RADIX == 16
 136   mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_16_DIGITS) * 4);
 137 #else
 138 #error "GFC_REAL_16_RADIX has unknown value"
 139 #endif
 140   v2 = v2 & mask;
 141   *f = (GFC_REAL_16) v1 * GFC_REAL_16_LITERAL(0x1.p-64)
 142     + (GFC_REAL_16) v2 * GFC_REAL_16_LITERAL(0x1.p-128);
 143 }
 144 #endif
 145 /* libgfortran previously had a Mersenne Twister, taken from the paper:
 146
 147         Mersenne Twister:       623-dimensionally equidistributed
 148                                 uniform pseudorandom generator.
 149
 150         by Makoto Matsumoto & Takuji Nishimura
 151         which appeared in the: ACM Transactions on Modelling and Computer
 152         Simulations: Special Issue on Uniform Random Number
 153         Generation. ( Early in 1998 ).
 154
 155    The Mersenne Twister code was replaced due to
 156
 157     (1) Simple user specified seeds lead to really bad sequences for
 158         nearly 100000 random numbers.
 159     (2) open(), read(), and close() were not properly declared via header
 160         files.
 161     (3) The global index i was abused and caused unexpected behavior with
 162         GET and PUT.
 163     (4) See PR 15619.
 164
 165
 166    libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
 167    random number generator.  This PRNG combines:
 168
 169    (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
 170        of 2^32,
 171    (2) A 3-shift shift-register generator with a period of 2^32-1,
 172    (3) Two 16-bit multiply-with-carry generators with a period of
 173        597273182964842497 > 2^59.
 174
 175    The overall period exceeds 2^123.
 176
 177    http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
 178
 179    The above web site has an archive of a newsgroup posting from George
 180    Marsaglia with the statement:
 181
 182    Subject: Random numbers for C: Improvements.
 183    Date: Fri, 15 Jan 1999 11:41:47 -0500
 184    From: George Marsaglia <geo@stat.fsu.edu>
 185    Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
 186    References: <369B5E30.65A55FD1@stat.fsu.edu>
 187    Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
 188    Lines: 93
 189
 190    As I hoped, several suggestions have led to
 191    improvements in the code for RNG's I proposed for
 192    use in C. (See the thread "Random numbers for C: Some
 193    suggestions" in previous postings.) The improved code
 194    is listed below.
 195
 196    A question of copyright has also been raised.  Unlike
 197    DIEHARD, there is no copyright on the code below. You
 198    are free to use it in any way you want, but you may
 199    wish to acknowledge the source, as a courtesy.
 200
 201 "There is no copyright on the code below." included the original
 202 KISS algorithm.  */
 203
 204 /* We use three KISS random number generators, with different
 205    seeds.
 206    As a matter of Quality of Implementation, the random numbers
 207    we generate for different REAL kinds, starting from the same
 208    seed, are always the same up to the precision of these types.
 209    We do this by using three generators with different seeds, the
 210    first one always for the most significant bits, the second one
 211    for bits 33..64 (if present in the REAL kind), and the third one
 212    (called twice) for REAL(16).  */
 213
 214 #define GFC_SL(k, n)    ((k)^((k)<<(n)))
 215 #define GFC_SR(k, n)    ((k)^((k)>>(n)))
 216
 217 /*   Reference for the seed:
 218    From: "George Marsaglia" <g...@stat.fsu.edu>
 219    Newsgroups: sci.math
 220    Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
 221
 222    The KISS RNG uses four seeds, x, y, z, c,
 223    with 0<=x<2^32, 0<y<2^32, 0<=z<2^32, 0<=c<698769069
 224    except that the two pairs
 225    z=0,c=0 and z=2^32-1,c=698769068
 226    should be avoided.  */
 227
 228 /* Any modifications to the seeds that change kiss_size below need to be
 229    reflected in check.c (gfc_check_random_seed) to enable correct
 230    compile-time checking of PUT size for the RANDOM_SEED intrinsic.  */
 231
 232 #define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069
 233 #define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021
 234 #ifdef HAVE_GFC_REAL_16
 235 #define KISS_DEFAULT_SEED_3 573658661, 185639104, 582619469, 296736107
 236 #endif
 237
 238 static GFC_UINTEGER_4 kiss_seed[] = {
 239   KISS_DEFAULT_SEED_1,
 240   KISS_DEFAULT_SEED_2,
 241 #ifdef HAVE_GFC_REAL_16
 242   KISS_DEFAULT_SEED_3
 243 #endif
 244 };
 245
 246 static GFC_UINTEGER_4 kiss_default_seed[] = {
 247   KISS_DEFAULT_SEED_1,
 248   KISS_DEFAULT_SEED_2,
 249 #ifdef HAVE_GFC_REAL_16
 250   KISS_DEFAULT_SEED_3
 251 #endif
 252 };
 253
 254 static const GFC_INTEGER_4 kiss_size = sizeof(kiss_seed)/sizeof(kiss_seed[0]);
 255
 256 static GFC_UINTEGER_4 * const kiss_seed_1 = kiss_seed;
 257 static GFC_UINTEGER_4 * const kiss_seed_2 = kiss_seed + 4;
 258
 259 #ifdef HAVE_GFC_REAL_16
 260 static GFC_UINTEGER_4 * const kiss_seed_3 = kiss_seed + 8;
 261 #endif
 262
 263 /* kiss_random_kernel() returns an integer value in the range of
 264    (0, GFC_UINTEGER_4_HUGE].  The distribution of pseudorandom numbers
 265    should be uniform.  */
 266
 267 static GFC_UINTEGER_4
 268 kiss_random_kernel(GFC_UINTEGER_4 * seed)
 269 {
 270   GFC_UINTEGER_4 kiss;
 271
 272   seed[0] = 69069 * seed[0] + 1327217885;
 273   seed[1] = GFC_SL(GFC_SR(GFC_SL(seed[1],13),17),5);
 274   seed[2] = 18000 * (seed[2] & 65535) + (seed[2] >> 16);
 275   seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16);
 276   kiss = seed[0] + seed[1] + (seed[2] << 16) + seed[3];
 277
 278   return kiss;
 279 }
 280
 281 /*  This function produces a REAL(4) value from the uniform distribution
 282     with range [0,1).  */
 283
 284 void
 285 random_r4 (GFC_REAL_4 *x)
 286 {
 287   GFC_UINTEGER_4 kiss;
 288
 289   __gthread_mutex_lock (&random_lock);
 290   kiss = kiss_random_kernel (kiss_seed_1);
 291   rnumber_4 (x, kiss);
 292   __gthread_mutex_unlock (&random_lock);
 293 }
 294 iexport(random_r4);
 295
 296 /*  This function produces a REAL(8) value from the uniform distribution
 297     with range [0,1).  */
 298
 299 void
 300 random_r8 (GFC_REAL_8 *x)
 301 {
 302   GFC_UINTEGER_8 kiss;
 303
 304   __gthread_mutex_lock (&random_lock);
 305   kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
 306   kiss += kiss_random_kernel (kiss_seed_2);
 307   rnumber_8 (x, kiss);
 308   __gthread_mutex_unlock (&random_lock);
 309 }
 310 iexport(random_r8);
 311
 312 #ifdef HAVE_GFC_REAL_10
 313
 314 /*  This function produces a REAL(10) value from the uniform distribution
 315     with range [0,1).  */
 316
 317 void
 318 random_r10 (GFC_REAL_10 *x)
 319 {
 320   GFC_UINTEGER_8 kiss;
 321
 322   __gthread_mutex_lock (&random_lock);
 323   kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
 324   kiss += kiss_random_kernel (kiss_seed_2);
 325   rnumber_10 (x, kiss);
 326   __gthread_mutex_unlock (&random_lock);
 327 }
 328 iexport(random_r10);
 329
 330 #endif
 331
 332 /*  This function produces a REAL(16) value from the uniform distribution
 333     with range [0,1).  */
 334
 335 #ifdef HAVE_GFC_REAL_16
 336
 337 void
 338 random_r16 (GFC_REAL_16 *x)
 339 {
 340   GFC_UINTEGER_8 kiss1, kiss2;
 341
 342   __gthread_mutex_lock (&random_lock);
 343   kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
 344   kiss1 += kiss_random_kernel (kiss_seed_2);
 345
 346   kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
 347   kiss2 += kiss_random_kernel (kiss_seed_3);
 348
 349   rnumber_16 (x, kiss1, kiss2);
 350   __gthread_mutex_unlock (&random_lock);
 351 }
 352 iexport(random_r16);
 353
 354
 355 #endif
 356 /*  This function fills a REAL(4) array with values from the uniform
 357     distribution with range [0,1).  */
 358
 359 void
 360 arandom_r4 (gfc_array_r4 *x)
 361 {
 362   index_type count[GFC_MAX_DIMENSIONS];
 363   index_type extent[GFC_MAX_DIMENSIONS];
 364   index_type stride[GFC_MAX_DIMENSIONS];
 365   index_type stride0;
 366   index_type dim;
 367   GFC_REAL_4 *dest;
 368   GFC_UINTEGER_4 kiss;
 369   int n;
 370
 371   dest = x->data;
 372
 373   dim = GFC_DESCRIPTOR_RANK (x);
 374
 375   for (n = 0; n < dim; n++)
 376     {
 377       count[n] = 0;
 378       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 379       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 380       if (extent[n] <= 0)
 381         return;
 382     }
 383
 384   stride0 = stride[0];
 385
 386   __gthread_mutex_lock (&random_lock);
 387
 388   while (dest)
 389     {
 390       /* random_r4 (dest);  */
 391       kiss = kiss_random_kernel (kiss_seed_1);
 392       rnumber_4 (dest, kiss);
 393
 394       /* Advance to the next element.  */
 395       dest += stride0;
 396       count[0]++;
 397       /* Advance to the next source element.  */
 398       n = 0;
 399       while (count[n] == extent[n])
 400         {
 401           /* When we get to the end of a dimension, reset it and increment
 402              the next dimension.  */
 403           count[n] = 0;
 404           /* We could precalculate these products, but this is a less
 405              frequently used path so probably not worth it.  */
 406           dest -= stride[n] * extent[n];
 407           n++;
 408           if (n == dim)
 409             {
 410               dest = NULL;
 411               break;
 412             }
 413           else
 414             {
 415               count[n]++;
 416               dest += stride[n];
 417             }
 418         }
 419     }
 420   __gthread_mutex_unlock (&random_lock);
 421 }
 422
 423 /*  This function fills a REAL(8) array with values from the uniform
 424     distribution with range [0,1).  */
 425
 426 void
 427 arandom_r8 (gfc_array_r8 *x)
 428 {
 429   index_type count[GFC_MAX_DIMENSIONS];
 430   index_type extent[GFC_MAX_DIMENSIONS];
 431   index_type stride[GFC_MAX_DIMENSIONS];
 432   index_type stride0;
 433   index_type dim;
 434   GFC_REAL_8 *dest;
 435   GFC_UINTEGER_8 kiss;
 436   int n;
 437
 438   dest = x->data;
 439
 440   dim = GFC_DESCRIPTOR_RANK (x);
 441
 442   for (n = 0; n < dim; n++)
 443     {
 444       count[n] = 0;
 445       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 446       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 447       if (extent[n] <= 0)
 448         return;
 449     }
 450
 451   stride0 = stride[0];
 452
 453   __gthread_mutex_lock (&random_lock);
 454
 455   while (dest)
 456     {
 457       /* random_r8 (dest);  */
 458       kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
 459       kiss += kiss_random_kernel (kiss_seed_2);
 460       rnumber_8 (dest, kiss);
 461
 462       /* Advance to the next element.  */
 463       dest += stride0;
 464       count[0]++;
 465       /* Advance to the next source element.  */
 466       n = 0;
 467       while (count[n] == extent[n])
 468         {
 469           /* When we get to the end of a dimension, reset it and increment
 470              the next dimension.  */
 471           count[n] = 0;
 472           /* We could precalculate these products, but this is a less
 473              frequently used path so probably not worth it.  */
 474           dest -= stride[n] * extent[n];
 475           n++;
 476           if (n == dim)
 477             {
 478               dest = NULL;
 479               break;
 480             }
 481           else
 482             {
 483               count[n]++;
 484               dest += stride[n];
 485             }
 486         }
 487     }
 488   __gthread_mutex_unlock (&random_lock);
 489 }
 490
 491 #ifdef HAVE_GFC_REAL_10
 492
 493 /*  This function fills a REAL(10) array with values from the uniform
 494     distribution with range [0,1).  */
 495
 496 void
 497 arandom_r10 (gfc_array_r10 *x)
 498 {
 499   index_type count[GFC_MAX_DIMENSIONS];
 500   index_type extent[GFC_MAX_DIMENSIONS];
 501   index_type stride[GFC_MAX_DIMENSIONS];
 502   index_type stride0;
 503   index_type dim;
 504   GFC_REAL_10 *dest;
 505   GFC_UINTEGER_8 kiss;
 506   int n;
 507
 508   dest = x->data;
 509
 510   dim = GFC_DESCRIPTOR_RANK (x);
 511
 512   for (n = 0; n < dim; n++)
 513     {
 514       count[n] = 0;
 515       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 516       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 517       if (extent[n] <= 0)
 518         return;
 519     }
 520
 521   stride0 = stride[0];
 522
 523   __gthread_mutex_lock (&random_lock);
 524
 525   while (dest)
 526     {
 527       /* random_r10 (dest);  */
 528       kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
 529       kiss += kiss_random_kernel (kiss_seed_2);
 530       rnumber_10 (dest, kiss);
 531
 532       /* Advance to the next element.  */
 533       dest += stride0;
 534       count[0]++;
 535       /* Advance to the next source element.  */
 536       n = 0;
 537       while (count[n] == extent[n])
 538         {
 539           /* When we get to the end of a dimension, reset it and increment
 540              the next dimension.  */
 541           count[n] = 0;
 542           /* We could precalculate these products, but this is a less
 543              frequently used path so probably not worth it.  */
 544           dest -= stride[n] * extent[n];
 545           n++;
 546           if (n == dim)
 547             {
 548               dest = NULL;
 549               break;
 550             }
 551           else
 552             {
 553               count[n]++;
 554               dest += stride[n];
 555             }
 556         }
 557     }
 558   __gthread_mutex_unlock (&random_lock);
 559 }
 560
 561 #endif
 562
 563 #ifdef HAVE_GFC_REAL_16
 564
 565 /*  This function fills a REAL(16) array with values from the uniform
 566     distribution with range [0,1).  */
 567
 568 void
 569 arandom_r16 (gfc_array_r16 *x)
 570 {
 571   index_type count[GFC_MAX_DIMENSIONS];
 572   index_type extent[GFC_MAX_DIMENSIONS];
 573   index_type stride[GFC_MAX_DIMENSIONS];
 574   index_type stride0;
 575   index_type dim;
 576   GFC_REAL_16 *dest;
 577   GFC_UINTEGER_8 kiss1, kiss2;
 578   int n;
 579
 580   dest = x->data;
 581
 582   dim = GFC_DESCRIPTOR_RANK (x);
 583
 584   for (n = 0; n < dim; n++)
 585     {
 586       count[n] = 0;
 587       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 588       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 589       if (extent[n] <= 0)
 590         return;
 591     }
 592
 593   stride0 = stride[0];
 594
 595   __gthread_mutex_lock (&random_lock);
 596
 597   while (dest)
 598     {
 599       /* random_r16 (dest);  */
 600       kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
 601       kiss1 += kiss_random_kernel (kiss_seed_2);
 602
 603       kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
 604       kiss2 += kiss_random_kernel (kiss_seed_3);
 605
 606       rnumber_16 (dest, kiss1, kiss2);
 607
 608       /* Advance to the next element.  */
 609       dest += stride0;
 610       count[0]++;
 611       /* Advance to the next source element.  */
 612       n = 0;
 613       while (count[n] == extent[n])
 614         {
 615           /* When we get to the end of a dimension, reset it and increment
 616              the next dimension.  */
 617           count[n] = 0;
 618           /* We could precalculate these products, but this is a less
 619              frequently used path so probably not worth it.  */
 620           dest -= stride[n] * extent[n];
 621           n++;
 622           if (n == dim)
 623             {
 624               dest = NULL;
 625               break;
 626             }
 627           else
 628             {
 629               count[n]++;
 630               dest += stride[n];
 631             }
 632         }
 633     }
 634   __gthread_mutex_unlock (&random_lock);
 635 }
 636
 637 #endif
 638
 639
 640
 641 static void
 642 scramble_seed (unsigned char *dest, unsigned char *src, int size)
 643 {
 644   int i;
 645
 646   for (i = 0; i < size; i++)
 647     dest[(i % 2) * (size / 2) + i / 2] = src[i];
 648 }
 649
 650
 651 static void
 652 unscramble_seed (unsigned char *dest, unsigned char *src, int size)
 653 {
 654   int i;
 655
 656   for (i = 0; i < size; i++)
 657     dest[i] = src[(i % 2) * (size / 2) + i / 2];
 658 }
 659
 660
 661
 662 /* random_seed is used to seed the PRNG with either a default
 663    set of seeds or user specified set of seeds.  random_seed
 664    must be called with no argument or exactly one argument.  */
 665
 666 void
 667 random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 668 {
 669   int i;
 670   unsigned char seed[4*kiss_size];
 671
 672   __gthread_mutex_lock (&random_lock);
 673
 674   /* Check that we only have one argument present.  */
 675   if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
 676     runtime_error ("RANDOM_SEED should have at most one argument present.");
 677
 678   /* From the standard: "If no argument is present, the processor assigns
 679      a processor-dependent value to the seed."  */
 680   if (size == NULL && put == NULL && get == NULL)
 681       for (i = 0; i < kiss_size; i++)
 682         kiss_seed[i] = kiss_default_seed[i];
 683
 684   if (size != NULL)
 685     *size = kiss_size;
 686
 687   if (put != NULL)
 688     {
 689       /* If the rank of the array is not 1, abort.  */
 690       if (GFC_DESCRIPTOR_RANK (put) != 1)
 691         runtime_error ("Array rank of PUT is not 1.");
 692
 693       /* If the array is too small, abort.  */
 694       if (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size)
 695         runtime_error ("Array size of PUT is too small.");
 696
 697       /*  We copy the seed given by the user.  */
 698       for (i = 0; i < kiss_size; i++)
 699         memcpy (seed + i * sizeof(GFC_UINTEGER_4),
 700                 &(put->data[(kiss_size - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]),
 701                 sizeof(GFC_UINTEGER_4));
 702
 703       /* We put it after scrambling the bytes, to paper around users who
 704          provide seeds with quality only in the lower or upper part.  */
 705       scramble_seed ((unsigned char *) kiss_seed, seed, 4*kiss_size);
 706     }
 707
 708   /* Return the seed to GET data.  */
 709   if (get != NULL)
 710     {
 711       /* If the rank of the array is not 1, abort.  */
 712       if (GFC_DESCRIPTOR_RANK (get) != 1)
 713         runtime_error ("Array rank of GET is not 1.");
 714
 715       /* If the array is too small, abort.  */
 716       if (GFC_DESCRIPTOR_EXTENT(get,0) < kiss_size)
 717         runtime_error ("Array size of GET is too small.");
 718
 719       /* Unscramble the seed.  */
 720       unscramble_seed (seed, (unsigned char *) kiss_seed, 4*kiss_size);
 721
 722       /*  Then copy it back to the user variable.  */
 723       for (i = 0; i < kiss_size; i++)
 724         memcpy (&(get->data[(kiss_size - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]),
 725                seed + i * sizeof(GFC_UINTEGER_4),
 726                sizeof(GFC_UINTEGER_4));
 727     }
 728
 729   __gthread_mutex_unlock (&random_lock);
 730 }
 731 iexport(random_seed_i4);
 732
 733
 734 void
 735 random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
 736 {
 737   int i;
 738
 739   __gthread_mutex_lock (&random_lock);
 740
 741   /* Check that we only have one argument present.  */
 742   if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
 743     runtime_error ("RANDOM_SEED should have at most one argument present.");
 744
 745   /* From the standard: "If no argument is present, the processor assigns
 746      a processor-dependent value to the seed."  */
 747   if (size == NULL && put == NULL && get == NULL)
 748       for (i = 0; i < kiss_size; i++)
 749         kiss_seed[i] = kiss_default_seed[i];
 750
 751   if (size != NULL)
 752     *size = kiss_size / 2;
 753
 754   if (put != NULL)
 755     {
 756       /* If the rank of the array is not 1, abort.  */
 757       if (GFC_DESCRIPTOR_RANK (put) != 1)
 758         runtime_error ("Array rank of PUT is not 1.");
 759
 760       /* If the array is too small, abort.  */
 761       if (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size / 2)
 762         runtime_error ("Array size of PUT is too small.");
 763
 764       /*  This code now should do correct strides.  */
 765       for (i = 0; i < kiss_size / 2; i++)
 766         memcpy (&kiss_seed[2*i], &(put->data[i * GFC_DESCRIPTOR_STRIDE(put,0)]),
 767                 sizeof (GFC_UINTEGER_8));
 768     }
 769
 770   /* Return the seed to GET data.  */
 771   if (get != NULL)
 772     {
 773       /* If the rank of the array is not 1, abort.  */
 774       if (GFC_DESCRIPTOR_RANK (get) != 1)
 775         runtime_error ("Array rank of GET is not 1.");
 776
 777       /* If the array is too small, abort.  */
 778       if (GFC_DESCRIPTOR_EXTENT(get,0) < kiss_size / 2)
 779         runtime_error ("Array size of GET is too small.");
 780
 781       /*  This code now should do correct strides.  */
 782       for (i = 0; i < kiss_size / 2; i++)
 783         memcpy (&(get->data[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &kiss_seed[2*i],
 784                 sizeof (GFC_UINTEGER_8));
 785     }
 786
 787   __gthread_mutex_unlock (&random_lock);
 788 }
 789 iexport(random_seed_i8);
 790
 791
 792 #ifndef __GTHREAD_MUTEX_INIT
 793 static void __attribute__((constructor))
 794 init (void)
 795 {
 796   __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
 797 }
 798 #endif