OSDN Git Service

cbe4f90322b8b37398efda06f27fec6dd743ac72
[pf3gnuchains/gcc-fork.git] / 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_LITERAL(0x1.p-32);
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_LITERAL(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_LITERAL(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_LITERAL(0x1.p-64)
141     + (GFC_REAL_16) v2 * GFC_REAL_16_LITERAL(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] = GFC_DESCRIPTOR_STRIDE(x,n);
378       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
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] = GFC_DESCRIPTOR_STRIDE(x,n);
445       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
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] = GFC_DESCRIPTOR_STRIDE(x,n);
515       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
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] = GFC_DESCRIPTOR_STRIDE(x,n);
587       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
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 (GFC_DESCRIPTOR_EXTENT(put,0) < 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) * GFC_DESCRIPTOR_STRIDE(put,0)]),
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 (GFC_DESCRIPTOR_EXTENT(get,0) < 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) * GFC_DESCRIPTOR_STRIDE(get,0)]),
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 (GFC_DESCRIPTOR_EXTENT(put,0) < 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 * GFC_DESCRIPTOR_STRIDE(put,0)]),
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 (GFC_DESCRIPTOR_EXTENT(get,0) < 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 * GFC_DESCRIPTOR_STRIDE(get,0)]), &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