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