libgfortran/intrinsics/random.c

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