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