* src/alpha/ffi.c: Do not include stdlib.h.
[official-gcc.git] / libgfortran / intrinsics / random.c
blob593a6518cae3b6b03111872b35283f36f37dbd0b
1 /* Implementation of the RANDOM intrinsics
2 Copyright (C) 2002-2014 Free Software Foundation, Inc.
3 Contributed by Lars Segerlund <seger@linuxmail.org>
4 and Steve Kargl.
6 This file is part of the GNU Fortran runtime library (libgfortran).
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.
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.
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.
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/>. */
27 #include "libgfortran.h"
28 #include <gthr.h>
29 #include <string.h>
31 extern void random_r4 (GFC_REAL_4 *);
32 iexport_proto(random_r4);
34 extern void random_r8 (GFC_REAL_8 *);
35 iexport_proto(random_r8);
37 extern void arandom_r4 (gfc_array_r4 *);
38 export_proto(arandom_r4);
40 extern void arandom_r8 (gfc_array_r8 *);
41 export_proto(arandom_r8);
43 #ifdef HAVE_GFC_REAL_10
45 extern void random_r10 (GFC_REAL_10 *);
46 iexport_proto(random_r10);
48 extern void arandom_r10 (gfc_array_r10 *);
49 export_proto(arandom_r10);
51 #endif
53 #ifdef HAVE_GFC_REAL_16
55 extern void random_r16 (GFC_REAL_16 *);
56 iexport_proto(random_r16);
58 extern void arandom_r16 (gfc_array_r16 *);
59 export_proto(arandom_r16);
61 #endif
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
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. */
76 static void
77 rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v)
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_LITERAL(0x1.p-32);
91 static void
92 rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v)
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_LITERAL(0x1.p-64);
106 #ifdef HAVE_GFC_REAL_10
108 static void
109 rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v)
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_LITERAL(0x1.p-64);
122 #endif
124 #ifdef HAVE_GFC_REAL_16
126 /* For REAL(KIND=16), we only need to mask off the lower bits. */
128 static void
129 rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
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_LITERAL(0x1.p-64)
141 + (GFC_REAL_16) v2 * GFC_REAL_16_LITERAL(0x1.p-128);
143 #endif
144 /* libgfortran previously had a Mersenne Twister, taken from the paper:
146 Mersenne Twister: 623-dimensionally equidistributed
147 uniform pseudorandom generator.
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 ).
154 The Mersenne Twister code was replaced due to
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.
165 libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
166 random number generator. This PRNG combines:
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.
174 The overall period exceeds 2^123.
176 http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
178 The above web site has an archive of a newsgroup posting from George
179 Marsaglia with the statement:
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
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.
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.
200 "There is no copyright on the code below." included the original
201 KISS algorithm. */
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). */
213 #define GFC_SL(k, n) ((k)^((k)<<(n)))
214 #define GFC_SR(k, n) ((k)^((k)>>(n)))
216 /* Reference for the seed:
217 From: "George Marsaglia" <g...@stat.fsu.edu>
218 Newsgroups: sci.math
219 Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
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. */
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. */
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
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
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
253 static const GFC_INTEGER_4 kiss_size = sizeof(kiss_seed)/sizeof(kiss_seed[0]);
255 static GFC_UINTEGER_4 * const kiss_seed_1 = kiss_seed;
256 static GFC_UINTEGER_4 * const kiss_seed_2 = kiss_seed + 4;
258 #ifdef HAVE_GFC_REAL_16
259 static GFC_UINTEGER_4 * const kiss_seed_3 = kiss_seed + 8;
260 #endif
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. */
266 static GFC_UINTEGER_4
267 kiss_random_kernel(GFC_UINTEGER_4 * seed)
269 GFC_UINTEGER_4 kiss;
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];
277 return kiss;
280 /* This function produces a REAL(4) value from the uniform distribution
281 with range [0,1). */
283 void
284 random_r4 (GFC_REAL_4 *x)
286 GFC_UINTEGER_4 kiss;
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);
293 iexport(random_r4);
295 /* This function produces a REAL(8) value from the uniform distribution
296 with range [0,1). */
298 void
299 random_r8 (GFC_REAL_8 *x)
301 GFC_UINTEGER_8 kiss;
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);
309 iexport(random_r8);
311 #ifdef HAVE_GFC_REAL_10
313 /* This function produces a REAL(10) value from the uniform distribution
314 with range [0,1). */
316 void
317 random_r10 (GFC_REAL_10 *x)
319 GFC_UINTEGER_8 kiss;
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);
327 iexport(random_r10);
329 #endif
331 /* This function produces a REAL(16) value from the uniform distribution
332 with range [0,1). */
334 #ifdef HAVE_GFC_REAL_16
336 void
337 random_r16 (GFC_REAL_16 *x)
339 GFC_UINTEGER_8 kiss1, kiss2;
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);
345 kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
346 kiss2 += kiss_random_kernel (kiss_seed_3);
348 rnumber_16 (x, kiss1, kiss2);
349 __gthread_mutex_unlock (&random_lock);
351 iexport(random_r16);
354 #endif
355 /* This function fills a REAL(4) array with values from the uniform
356 distribution with range [0,1). */
358 void
359 arandom_r4 (gfc_array_r4 *x)
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;
370 dest = x->base_addr;
372 dim = GFC_DESCRIPTOR_RANK (x);
374 for (n = 0; n < dim; n++)
376 count[n] = 0;
377 stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
378 extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
379 if (extent[n] <= 0)
380 return;
383 stride0 = stride[0];
385 __gthread_mutex_lock (&random_lock);
387 while (dest)
389 /* random_r4 (dest); */
390 kiss = kiss_random_kernel (kiss_seed_1);
391 rnumber_4 (dest, kiss);
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])
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)
409 dest = NULL;
410 break;
412 else
414 count[n]++;
415 dest += stride[n];
419 __gthread_mutex_unlock (&random_lock);
422 /* This function fills a REAL(8) array with values from the uniform
423 distribution with range [0,1). */
425 void
426 arandom_r8 (gfc_array_r8 *x)
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;
437 dest = x->base_addr;
439 dim = GFC_DESCRIPTOR_RANK (x);
441 for (n = 0; n < dim; n++)
443 count[n] = 0;
444 stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
445 extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
446 if (extent[n] <= 0)
447 return;
450 stride0 = stride[0];
452 __gthread_mutex_lock (&random_lock);
454 while (dest)
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);
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])
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)
477 dest = NULL;
478 break;
480 else
482 count[n]++;
483 dest += stride[n];
487 __gthread_mutex_unlock (&random_lock);
490 #ifdef HAVE_GFC_REAL_10
492 /* This function fills a REAL(10) array with values from the uniform
493 distribution with range [0,1). */
495 void
496 arandom_r10 (gfc_array_r10 *x)
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;
507 dest = x->base_addr;
509 dim = GFC_DESCRIPTOR_RANK (x);
511 for (n = 0; n < dim; n++)
513 count[n] = 0;
514 stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
515 extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
516 if (extent[n] <= 0)
517 return;
520 stride0 = stride[0];
522 __gthread_mutex_lock (&random_lock);
524 while (dest)
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);
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])
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)
547 dest = NULL;
548 break;
550 else
552 count[n]++;
553 dest += stride[n];
557 __gthread_mutex_unlock (&random_lock);
560 #endif
562 #ifdef HAVE_GFC_REAL_16
564 /* This function fills a REAL(16) array with values from the uniform
565 distribution with range [0,1). */
567 void
568 arandom_r16 (gfc_array_r16 *x)
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;
579 dest = x->base_addr;
581 dim = GFC_DESCRIPTOR_RANK (x);
583 for (n = 0; n < dim; n++)
585 count[n] = 0;
586 stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
587 extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
588 if (extent[n] <= 0)
589 return;
592 stride0 = stride[0];
594 __gthread_mutex_lock (&random_lock);
596 while (dest)
598 /* random_r16 (dest); */
599 kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
600 kiss1 += kiss_random_kernel (kiss_seed_2);
602 kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
603 kiss2 += kiss_random_kernel (kiss_seed_3);
605 rnumber_16 (dest, kiss1, kiss2);
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])
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)
623 dest = NULL;
624 break;
626 else
628 count[n]++;
629 dest += stride[n];
633 __gthread_mutex_unlock (&random_lock);
636 #endif
640 static void
641 scramble_seed (unsigned char *dest, unsigned char *src, int size)
643 int i;
645 for (i = 0; i < size; i++)
646 dest[(i % 2) * (size / 2) + i / 2] = src[i];
650 static void
651 unscramble_seed (unsigned char *dest, unsigned char *src, int size)
653 int i;
655 for (i = 0; i < size; i++)
656 dest[i] = src[(i % 2) * (size / 2) + i / 2];
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. */
665 void
666 random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
668 int i;
669 unsigned char seed[4*kiss_size];
671 __gthread_mutex_lock (&random_lock);
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.");
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];
683 if (size != NULL)
684 *size = kiss_size;
686 if (put != NULL)
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.");
692 /* If the array is too small, abort. */
693 if (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size)
694 runtime_error ("Array size of PUT is too small.");
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->base_addr[(kiss_size - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]),
700 sizeof(GFC_UINTEGER_4));
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);
707 /* Return the seed to GET data. */
708 if (get != NULL)
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.");
714 /* If the array is too small, abort. */
715 if (GFC_DESCRIPTOR_EXTENT(get,0) < kiss_size)
716 runtime_error ("Array size of GET is too small.");
718 /* Unscramble the seed. */
719 unscramble_seed (seed, (unsigned char *) kiss_seed, 4*kiss_size);
721 /* Then copy it back to the user variable. */
722 for (i = 0; i < kiss_size; i++)
723 memcpy (&(get->base_addr[(kiss_size - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]),
724 seed + i * sizeof(GFC_UINTEGER_4),
725 sizeof(GFC_UINTEGER_4));
728 __gthread_mutex_unlock (&random_lock);
730 iexport(random_seed_i4);
733 void
734 random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
736 int i;
738 __gthread_mutex_lock (&random_lock);
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.");
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];
750 if (size != NULL)
751 *size = kiss_size / 2;
753 if (put != NULL)
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.");
759 /* If the array is too small, abort. */
760 if (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size / 2)
761 runtime_error ("Array size of PUT is too small.");
763 /* This code now should do correct strides. */
764 for (i = 0; i < kiss_size / 2; i++)
765 memcpy (&kiss_seed[2*i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]),
766 sizeof (GFC_UINTEGER_8));
769 /* Return the seed to GET data. */
770 if (get != NULL)
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.");
776 /* If the array is too small, abort. */
777 if (GFC_DESCRIPTOR_EXTENT(get,0) < kiss_size / 2)
778 runtime_error ("Array size of GET is too small.");
780 /* This code now should do correct strides. */
781 for (i = 0; i < kiss_size / 2; i++)
782 memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &kiss_seed[2*i],
783 sizeof (GFC_UINTEGER_8));
786 __gthread_mutex_unlock (&random_lock);
788 iexport(random_seed_i8);
791 #ifndef __GTHREAD_MUTEX_INIT
792 static void __attribute__((constructor))
793 init (void)
795 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
797 #endif