OSDN Git Service

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