libgfortran/intrinsics/random.c

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