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 "config.h"
  33 #include "libgfortran.h"
  34 #include <gthr.h>
  35 #include <string.h>
  36
  37 extern void random_r4 (GFC_REAL_4 *);
  38 iexport_proto(random_r4);
  39
  40 extern void random_r8 (GFC_REAL_8 *);
  41 iexport_proto(random_r8);
  42
  43 extern void arandom_r4 (gfc_array_r4 *);
  44 export_proto(arandom_r4);
  45
  46 extern void arandom_r8 (gfc_array_r8 *);
  47 export_proto(arandom_r8);
  48
  49 #ifdef HAVE_GFC_REAL_10
  50
  51 extern void random_r10 (GFC_REAL_10 *);
  52 iexport_proto(random_r10);
  53
  54 extern void arandom_r10 (gfc_array_r10 *);
  55 export_proto(arandom_r10);
  56
  57 #endif
  58
  59 #ifdef HAVE_GFC_REAL_16
  60
  61 extern void random_r16 (GFC_REAL_16 *);
  62 iexport_proto(random_r16);
  63
  64 extern void arandom_r16 (gfc_array_r16 *);
  65 export_proto(arandom_r16);
  66
  67 #endif
  68
  69 #ifdef __GTHREAD_MUTEX_INIT
  70 static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
  71 #else
  72 static __gthread_mutex_t random_lock;
  73 #endif
  74
  75 /* Helper routines to map a GFC_UINTEGER_* to the corresponding
  76    GFC_REAL_* types in the range of [0,1).  If GFC_REAL_*_RADIX are 2
  77    or 16, respectively, we mask off the bits that don't fit into the
  78    correct GFC_REAL_*, convert to the real type, then multiply by the
  79    correct offset.
  80 */
  81
  82
  83 static inline void
  84 rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v)
  85 {
  86   GFC_UINTEGER_4 mask;
  87 #if GFC_REAL_4_RADIX == 2
  88   mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS);
  89 #elif GFC_REAL_4_RADIX == 16
  90   mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4);
  91 #else
  92 #error "GFC_REAL_4_RADIX has unknown value"
  93 #endif
  94   v = v & mask;
  95   *f = (GFC_REAL_4) v * (GFC_REAL_4) 0x1.p-32f;
  96 }
  97
  98 static inline void
  99 rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v)
 100 {
 101   GFC_UINTEGER_8 mask;
 102 #if GFC_REAL_8_RADIX == 2
 103   mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_8_DIGITS);
 104 #elif GFC_REAL_8_RADIX == 16
 105   mask = ~ (GFC_UINTEGER_8) 0u << (16 - GFC_REAL_8_DIGITS) * 4);
 106 #else
 107 #error "GFC_REAL_8_RADIX has unknown value"
 108 #endif
 109   v = v & mask;
 110   *f = (GFC_REAL_8) v * (GFC_REAL_8) 0x1.p-64;
 111 }
 112
 113 #ifdef HAVE_GFC_REAL_10
 114
 115 static inline void
 116 rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v)
 117 {
 118   GFC_UINTEGER_8 mask;
 119 #if GFC_REAL_10_RADIX == 2
 120   mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_10_DIGITS);
 121 #elif GFC_REAL_10_RADIX == 16
 122   mask = ~ (GFC_UINTEGER_10) 0u << ((16 - GFC_REAL_10_DIGITS) * 4);
 123 #else
 124 #error "GFC_REAL_10_RADIX has unknown value"
 125 #endif
 126   v = v & mask;
 127   *f = (GFC_REAL_10) v * (GFC_REAL_10) 0x1.p-64;
 128 }
 129 #endif
 130
 131 #ifdef HAVE_GFC_REAL_16
 132
 133 /* For REAL(KIND=16), we only need to mask off the lower bits.  */
 134
 135 static inline void
 136 rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
 137 {
 138   GFC_UINTEGER_8 mask;
 139 #if GFC_REAL_16_RADIX == 2
 140   mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_16_DIGITS);
 141 #elif GFC_REAL_16_RADIX == 16
 142   mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_16_DIGITS) * 4);
 143 #else
 144 #error "GFC_REAL_16_RADIX has unknown value"
 145 #endif
 146   v2 = v2 & mask;
 147   *f = (GFC_REAL_16) v1 * (GFC_REAL_16) 0x1.p-64
 148     + (GFC_REAL_16) v2 * (GFC_REAL_16) 0x1.p-128;
 149 }
 150 #endif
 151 /* libgfortran previously had a Mersenne Twister, taken from the paper:
 152
 153         Mersenne Twister:       623-dimensionally equidistributed
 154                                 uniform pseudorandom generator.
 155
 156         by Makoto Matsumoto & Takuji Nishimura
 157         which appeared in the: ACM Transactions on Modelling and Computer
 158         Simulations: Special Issue on Uniform Random Number
 159         Generation. ( Early in 1998 ).
 160
 161    The Mersenne Twister code was replaced due to
 162
 163     (1) Simple user specified seeds lead to really bad sequences for
 164         nearly 100000 random numbers.
 165     (2) open(), read(), and close() were not properly declared via header
 166         files.
 167     (3) The global index i was abused and caused unexpected behavior with
 168         GET and PUT.
 169     (4) See PR 15619.
 170
 171
 172    libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
 173    random number generator.  This PRNG combines:
 174
 175    (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
 176        of 2^32,
 177    (2) A 3-shift shift-register generator with a period of 2^32-1,
 178    (3) Two 16-bit multiply-with-carry generators with a period of
 179        597273182964842497 > 2^59.
 180
 181    The overall period exceeds 2^123.
 182
 183    http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
 184
 185    The above web site has an archive of a newsgroup posting from George
 186    Marsaglia with the statement:
 187
 188    Subject: Random numbers for C: Improvements.
 189    Date: Fri, 15 Jan 1999 11:41:47 -0500
 190    From: George Marsaglia <geo@stat.fsu.edu>
 191    Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
 192    References: <369B5E30.65A55FD1@stat.fsu.edu>
 193    Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
 194    Lines: 93
 195
 196    As I hoped, several suggestions have led to
 197    improvements in the code for RNG's I proposed for
 198    use in C. (See the thread "Random numbers for C: Some
 199    suggestions" in previous postings.) The improved code
 200    is listed below.
 201
 202    A question of copyright has also been raised.  Unlike
 203    DIEHARD, there is no copyright on the code below. You
 204    are free to use it in any way you want, but you may
 205    wish to acknowledge the source, as a courtesy.
 206
 207 "There is no copyright on the code below." included the original
 208 KISS algorithm.  */
 209
 210 /* We use three KISS random number generators, with different
 211    seeds.
 212    As a matter of Quality of Implementation, the random numbers
 213    we generate for different REAL kinds, starting from the same
 214    seed, are always the same up to the precision of these types.
 215    We do this by using three generators with different seeds, the
 216    first one always for the most significant bits, the second one
 217    for bits 33..64 (if present in the REAL kind), and the third one
 218    (called twice) for REAL(16).
 219 */
 220
 221 #define GFC_SL(k, n)    ((k)^((k)<<(n)))
 222 #define GFC_SR(k, n)    ((k)^((k)>>(n)))
 223
 224 /*   Reference for the seed:
 225    From: "George Marsaglia" <g...@stat.fsu.edu>
 226    Newsgroups: sci.math
 227    Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
 228
 229    The KISS RNG uses four seeds, x, y, z, c,
 230    with 0<=x<2^32, 0<y<2^32, 0<=z<2^32, 0<=c<698769069
 231    except that the two pairs
 232    z=0,c=0 and z=2^32-1,c=698769068
 233    should be avoided.
 234 */
 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 /* random_seed is used to seed the PRNG with either a default
 644    set of seeds or user specified set of seeds.  random_seed
 645    must be called with no argument or exactly one argument.  */
 646
 647 void
 648 random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 649 {
 650   int i;
 651
 652   __gthread_mutex_lock (&random_lock);
 653
 654   /* Check that we only have one argument present.  */
 655   if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
 656     runtime_error ("RANDOM_SEED should have at most one argument present.");
 657
 658   /* From the standard: "If no argument is present, the processor assigns
 659      a processor-dependent value to the seed."  */
 660   if (size == NULL && put == NULL && get == NULL)
 661       for (i = 0; i < kiss_size; i++)
 662         kiss_seed[i] = kiss_default_seed[i];
 663
 664   if (size != NULL)
 665     *size = kiss_size;
 666
 667   if (put != NULL)
 668     {
 669       /* If the rank of the array is not 1, abort.  */
 670       if (GFC_DESCRIPTOR_RANK (put) != 1)
 671         runtime_error ("Array rank of PUT is not 1.");
 672
 673       /* If the array is too small, abort.  */
 674       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
 675         runtime_error ("Array size of PUT is too small.");
 676
 677       /*  This code now should do correct strides.  */
 678       for (i = 0; i < kiss_size; i++)
 679         kiss_seed[i] = (GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
 680     }
 681
 682   /* Return the seed to GET data.  */
 683   if (get != NULL)
 684     {
 685       /* If the rank of the array is not 1, abort.  */
 686       if (GFC_DESCRIPTOR_RANK (get) != 1)
 687         runtime_error ("Array rank of GET is not 1.");
 688
 689       /* If the array is too small, abort.  */
 690       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
 691         runtime_error ("Array size of GET is too small.");
 692
 693       /*  This code now should do correct strides.  */
 694       for (i = 0; i < kiss_size; i++)
 695         get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
 696     }
 697
 698   __gthread_mutex_unlock (&random_lock);
 699 }
 700 iexport(random_seed_i4);
 701
 702
 703 void
 704 random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
 705 {
 706   int i;
 707
 708   __gthread_mutex_lock (&random_lock);
 709
 710   /* Check that we only have one argument present.  */
 711   if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
 712     runtime_error ("RANDOM_SEED should have at most one argument present.");
 713
 714   /* From the standard: "If no argument is present, the processor assigns
 715      a processor-dependent value to the seed."  */
 716   if (size == NULL && put == NULL && get == NULL)
 717       for (i = 0; i < kiss_size; i++)
 718         kiss_seed[i] = kiss_default_seed[i];
 719
 720   if (size != NULL)
 721     *size = kiss_size / 2;
 722
 723   if (put != NULL)
 724     {
 725       /* If the rank of the array is not 1, abort.  */
 726       if (GFC_DESCRIPTOR_RANK (put) != 1)
 727         runtime_error ("Array rank of PUT is not 1.");
 728
 729       /* If the array is too small, abort.  */
 730       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size / 2)
 731         runtime_error ("Array size of PUT is too small.");
 732
 733       /*  This code now should do correct strides.  */
 734       for (i = 0; i < kiss_size / 2; i++)
 735         memcpy (&kiss_seed[2*i], &(put->data[i * put->dim[0].stride]),
 736                 sizeof (GFC_UINTEGER_8));
 737     }
 738
 739   /* Return the seed to GET data.  */
 740   if (get != NULL)
 741     {
 742       /* If the rank of the array is not 1, abort.  */
 743       if (GFC_DESCRIPTOR_RANK (get) != 1)
 744         runtime_error ("Array rank of GET is not 1.");
 745
 746       /* If the array is too small, abort.  */
 747       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size / 2)
 748         runtime_error ("Array size of GET is too small.");
 749
 750       /*  This code now should do correct strides.  */
 751       for (i = 0; i < kiss_size / 2; i++)
 752         memcpy (&(get->data[i * get->dim[0].stride]), &kiss_seed[2*i],
 753                 sizeof (GFC_UINTEGER_8));
 754     }
 755
 756   __gthread_mutex_unlock (&random_lock);
 757 }
 758 iexport(random_seed_i8);
 759
 760
 761 #ifndef __GTHREAD_MUTEX_INIT
 762 static void __attribute__((constructor))
 763 init (void)
 764 {
 765   __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
 766 }
 767 #endif