OSDN Git Service

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