]>
Commit | Line | Data |
---|---|---|
6de9cd9a | 1 | /* Implementation of the RANDOM intrinsics |
4b6903ec | 2 | Copyright 2002, 2004, 2005 Free Software Foundation, Inc. |
6de9cd9a | 3 | Contributed by Lars Segerlund <seger@linuxmail.org> |
5f251c26 | 4 | and Steve Kargl. |
6de9cd9a DN |
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 | |
57dea9f6 | 9 | modify it under the terms of the GNU General Public |
6de9cd9a | 10 | License as published by the Free Software Foundation; either |
57dea9f6 TM |
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.) | |
6de9cd9a DN |
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 | |
57dea9f6 | 25 | GNU General Public License for more details. |
6de9cd9a | 26 | |
57dea9f6 TM |
27 | You should have received a copy of the GNU General Public |
28 | License along with libgfortran; see the file COPYING. If not, | |
fe2ae685 KC |
29 | write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, |
30 | Boston, MA 02110-1301, USA. */ | |
6de9cd9a | 31 | |
7d7b8bfe RH |
32 | #include "libgfortran.h" |
33 | ||
34 | extern void random_r4 (GFC_REAL_4 *); | |
35 | iexport_proto(random_r4); | |
36 | ||
37 | extern void random_r8 (GFC_REAL_8 *); | |
38 | iexport_proto(random_r8); | |
39 | ||
40 | extern void arandom_r4 (gfc_array_r4 *); | |
41 | export_proto(arandom_r4); | |
42 | ||
43 | extern void arandom_r8 (gfc_array_r8 *); | |
44 | export_proto(arandom_r8); | |
45 | ||
5f251c26 SK |
46 | #if 0 |
47 | ||
48 | /* The Mersenne Twister code is currently commented out due to | |
49 | ||
50 | (1) Simple user specified seeds lead to really bad sequences for | |
51 | nearly 100000 random numbers. | |
52 | (2) open(), read(), and close() are not properly declared via header | |
53 | files. | |
54 | (3) The global index i is abused and causes unexpected behavior with | |
55 | GET and PUT. | |
56 | (4) See PR 15619. | |
57 | ||
58 | The algorithm was taken from the paper : | |
59 | ||
60 | Mersenne Twister: 623-dimensionally equidistributed | |
61 | uniform pseudorandom generator. | |
62 | ||
63 | by: Makoto Matsumoto | |
64 | Takuji Nishimura | |
65 | ||
66 | Which appeared in the: ACM Transactions on Modelling and Computer | |
67 | Simulations: Special Issue on Uniform Random Number | |
68 | Generation. ( Early in 1998 ). */ | |
69 | ||
70 | ||
6de9cd9a DN |
71 | #include <stdio.h> |
72 | #include <stdlib.h> | |
73 | #include <sys/types.h> | |
74 | #include <sys/stat.h> | |
75 | #include <fcntl.h> | |
ebeb17c7 AJ |
76 | |
77 | #ifdef HAVE_UNISTD_H | |
78 | #include <unistd.h> | |
79 | #endif | |
80 | ||
6de9cd9a DN |
81 | /*Use the 'big' generator by default ( period -> 2**19937 ). */ |
82 | ||
83 | #define MT19937 | |
84 | ||
85 | /* Define the necessary constants for the algorithm. */ | |
86 | ||
87 | #ifdef MT19937 | |
88 | enum constants | |
89 | { | |
90 | N = 624, M = 397, R = 19, TU = 11, TS = 7, TT = 15, TL = 17 | |
91 | }; | |
92 | #define M_A 0x9908B0DF | |
93 | #define T_B 0x9D2C5680 | |
94 | #define T_C 0xEFC60000 | |
95 | #else | |
96 | enum constants | |
97 | { | |
98 | N = 351, M = 175, R = 19, TU = 11, TS = 7, TT = 15, TL = 17 | |
99 | }; | |
100 | #define M_A 0xE4BD75F5 | |
101 | #define T_B 0x655E5280 | |
102 | #define T_C 0xFFD58000 | |
103 | #endif | |
104 | ||
105 | static int i = N; | |
106 | static unsigned int seed[N]; | |
107 | ||
108 | /* This is the routine which handles the seeding of the generator, | |
109 | and also reading and writing of the seed. */ | |
110 | ||
6de9cd9a | 111 | void |
7d7b8bfe | 112 | random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get) |
6de9cd9a DN |
113 | { |
114 | /* Initialize the seed in system dependent manner. */ | |
115 | if (get == NULL && put == NULL && size == NULL) | |
116 | { | |
117 | int fd; | |
118 | fd = open ("/dev/urandom", O_RDONLY); | |
119 | if (fd == 0) | |
120 | { | |
121 | /* We dont have urandom. */ | |
122 | GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed; | |
123 | for (i = 0; i < N; i++) | |
124 | { | |
125 | s = s * 29943829 - 1; | |
126 | seed[i] = s; | |
127 | } | |
128 | } | |
129 | else | |
130 | { | |
131 | /* Using urandom, might have a length issue. */ | |
132 | read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N); | |
133 | close (fd); | |
134 | } | |
135 | return; | |
136 | } | |
137 | ||
138 | /* Return the size of the seed */ | |
139 | if (size != NULL) | |
140 | { | |
141 | *size = N; | |
142 | return; | |
143 | } | |
144 | ||
145 | /* if we have gotten to this pount we have a get or put | |
146 | * now we check it the array fulfills the demands in the standard . | |
147 | */ | |
148 | ||
149 | /* Set the seed to PUT data */ | |
150 | if (put != NULL) | |
151 | { | |
152 | /* if the rank of the array is not 1 abort */ | |
153 | if (GFC_DESCRIPTOR_RANK (put) != 1) | |
154 | abort (); | |
155 | ||
156 | /* if the array is too small abort */ | |
157 | if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N) | |
158 | abort (); | |
159 | ||
160 | /* If this is the case the array is a temporary */ | |
322c2425 | 161 | if (put->dim[0].stride == 0) |
6de9cd9a DN |
162 | return; |
163 | ||
164 | /* This code now should do correct strides. */ | |
165 | for (i = 0; i < N; i++) | |
166 | seed[i] = put->data[i * put->dim[0].stride]; | |
167 | } | |
168 | ||
169 | /* Return the seed to GET data */ | |
170 | if (get != NULL) | |
171 | { | |
172 | /* if the rank of the array is not 1 abort */ | |
173 | if (GFC_DESCRIPTOR_RANK (get) != 1) | |
174 | abort (); | |
175 | ||
176 | /* if the array is too small abort */ | |
177 | if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N) | |
178 | abort (); | |
179 | ||
180 | /* If this is the case the array is a temporary */ | |
181 | if (get->dim[0].stride == 0) | |
182 | return; | |
183 | ||
184 | /* This code now should do correct strides. */ | |
185 | for (i = 0; i < N; i++) | |
186 | get->data[i * get->dim[0].stride] = seed[i]; | |
187 | } | |
188 | } | |
7d7b8bfe | 189 | iexport(random_seed); |
6de9cd9a DN |
190 | |
191 | /* Here is the internal routine which generates the random numbers | |
192 | in 'batches' based upon the need for a new batch. | |
193 | It's an integer based routine known as 'Mersenne Twister'. | |
194 | This implementation still lacks 'tempering' and a good verification, | |
195 | but gives very good metrics. */ | |
196 | ||
197 | static void | |
198 | random_generate (void) | |
199 | { | |
200 | /* 32 bits. */ | |
201 | GFC_UINTEGER_4 y; | |
202 | ||
203 | /* Generate batch of N. */ | |
204 | int k, m; | |
205 | for (k = 0, m = M; k < N - 1; k++) | |
206 | { | |
207 | y = (seed[k] & (-1 << R)) | (seed[k + 1] & ((1u << R) - 1)); | |
208 | seed[k] = seed[m] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A); | |
209 | if (++m >= N) | |
210 | m = 0; | |
211 | } | |
212 | ||
213 | y = (seed[N - 1] & (-1 << R)) | (seed[0] & ((1u << R) - 1)); | |
214 | seed[N - 1] = seed[M - 1] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A); | |
215 | i = 0; | |
216 | } | |
217 | ||
218 | /* A routine to return a REAL(KIND=4). */ | |
219 | ||
6de9cd9a DN |
220 | void |
221 | random_r4 (GFC_REAL_4 * harv) | |
222 | { | |
223 | /* Regenerate if we need to. */ | |
224 | if (i >= N) | |
225 | random_generate (); | |
226 | ||
227 | /* Convert uint32 to REAL(KIND=4). */ | |
228 | *harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] / | |
229 | (GFC_REAL_4) (~(GFC_UINTEGER_4) 0)); | |
230 | } | |
7d7b8bfe | 231 | iexport(random_r4); |
6de9cd9a DN |
232 | |
233 | /* A routine to return a REAL(KIND=8). */ | |
234 | ||
6de9cd9a DN |
235 | void |
236 | random_r8 (GFC_REAL_8 * harv) | |
237 | { | |
238 | /* Regenerate if we need to, may waste one 32-bit value. */ | |
239 | if ((i + 1) >= N) | |
240 | random_generate (); | |
241 | ||
242 | /* Convert two uint32 to a REAL(KIND=8). */ | |
243 | *harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) / | |
244 | (GFC_REAL_8) (~(GFC_UINTEGER_8) 0); | |
245 | i += 2; | |
246 | } | |
7d7b8bfe | 247 | iexport(random_r8); |
6de9cd9a DN |
248 | |
249 | /* Code to handle arrays will follow here. */ | |
250 | ||
251 | /* REAL(KIND=4) REAL array. */ | |
252 | ||
6de9cd9a DN |
253 | void |
254 | arandom_r4 (gfc_array_r4 * harv) | |
255 | { | |
e33e218b TK |
256 | index_type count[GFC_MAX_DIMENSIONS]; |
257 | index_type extent[GFC_MAX_DIMENSIONS]; | |
258 | index_type stride[GFC_MAX_DIMENSIONS]; | |
6de9cd9a DN |
259 | index_type stride0; |
260 | index_type dim; | |
261 | GFC_REAL_4 *dest; | |
262 | int n; | |
263 | ||
264 | dest = harv->data; | |
265 | ||
266 | if (harv->dim[0].stride == 0) | |
267 | harv->dim[0].stride = 1; | |
268 | ||
269 | dim = GFC_DESCRIPTOR_RANK (harv); | |
270 | ||
271 | for (n = 0; n < dim; n++) | |
272 | { | |
273 | count[n] = 0; | |
274 | stride[n] = harv->dim[n].stride; | |
275 | extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound; | |
276 | if (extent[n] <= 0) | |
277 | return; | |
278 | } | |
279 | ||
280 | stride0 = stride[0]; | |
281 | ||
282 | while (dest) | |
283 | { | |
284 | /* Set the elements. */ | |
285 | ||
286 | /* regenerate if we need to */ | |
287 | if (i >= N) | |
288 | random_generate (); | |
289 | ||
290 | /* Convert uint32 to float in a hopefully g95 compiant manner */ | |
291 | *dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] / | |
292 | (GFC_REAL_4) (~(GFC_UINTEGER_4) 0)); | |
293 | ||
294 | /* Advance to the next element. */ | |
295 | dest += stride0; | |
296 | count[0]++; | |
297 | /* Advance to the next source element. */ | |
298 | n = 0; | |
299 | while (count[n] == extent[n]) | |
300 | { | |
301 | /* When we get to the end of a dimension, | |
302 | reset it and increment | |
303 | the next dimension. */ | |
304 | count[n] = 0; | |
305 | /* We could precalculate these products, | |
306 | but this is a less | |
307 | frequently used path so proabably not worth it. */ | |
308 | dest -= stride[n] * extent[n]; | |
309 | n++; | |
310 | if (n == dim) | |
311 | { | |
312 | dest = NULL; | |
313 | break; | |
314 | } | |
315 | else | |
316 | { | |
317 | count[n]++; | |
318 | dest += stride[n]; | |
319 | } | |
320 | } | |
321 | } | |
322 | } | |
323 | ||
324 | /* REAL(KIND=8) array. */ | |
325 | ||
6de9cd9a DN |
326 | void |
327 | arandom_r8 (gfc_array_r8 * harv) | |
328 | { | |
e33e218b TK |
329 | index_type count[GFC_MAX_DIMENSIONS]; |
330 | index_type extent[GFC_MAX_DIMENSIONS]; | |
331 | index_type stride[GFC_MAX_DIMENSIONS]; | |
6de9cd9a DN |
332 | index_type stride0; |
333 | index_type dim; | |
334 | GFC_REAL_8 *dest; | |
335 | int n; | |
336 | ||
337 | dest = harv->data; | |
338 | ||
339 | if (harv->dim[0].stride == 0) | |
340 | harv->dim[0].stride = 1; | |
341 | ||
342 | dim = GFC_DESCRIPTOR_RANK (harv); | |
343 | ||
344 | for (n = 0; n < dim; n++) | |
345 | { | |
346 | count[n] = 0; | |
347 | stride[n] = harv->dim[n].stride; | |
348 | extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound; | |
349 | if (extent[n] <= 0) | |
350 | return; | |
351 | } | |
352 | ||
353 | stride0 = stride[0]; | |
354 | ||
355 | while (dest) | |
356 | { | |
357 | /* Set the elements. */ | |
358 | ||
359 | /* regenerate if we need to, may waste one 32-bit value */ | |
360 | if ((i + 1) >= N) | |
361 | random_generate (); | |
362 | ||
363 | /* Convert two uint32 to a REAL(KIND=8). */ | |
364 | *dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) / | |
365 | (GFC_REAL_8) (~(GFC_UINTEGER_8) 0); | |
366 | i += 2; | |
367 | ||
368 | /* Advance to the next element. */ | |
369 | dest += stride0; | |
370 | count[0]++; | |
371 | /* Advance to the next source element. */ | |
372 | n = 0; | |
373 | while (count[n] == extent[n]) | |
374 | { | |
375 | /* When we get to the end of a dimension, | |
376 | reset it and increment | |
377 | the next dimension. */ | |
378 | count[n] = 0; | |
379 | /* We could precalculate these products, | |
380 | but this is a less | |
381 | frequently used path so proabably not worth it. */ | |
382 | dest -= stride[n] * extent[n]; | |
383 | n++; | |
384 | if (n == dim) | |
385 | { | |
386 | dest = NULL; | |
387 | break; | |
388 | } | |
389 | else | |
390 | { | |
391 | count[n]++; | |
392 | dest += stride[n]; | |
393 | } | |
394 | } | |
395 | } | |
396 | } | |
5f251c26 | 397 | |
7d7b8bfe | 398 | #else |
5f251c26 SK |
399 | |
400 | /* George Marsaglia's KISS (Keep It Simple Stupid) random number generator. | |
401 | ||
402 | This PRNG combines: | |
403 | ||
404 | (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period | |
405 | of 2^32, | |
406 | (2) A 3-shift shift-register generator with a period of 2^32-1, | |
407 | (3) Two 16-bit multiply-with-carry generators with a period of | |
408 | 597273182964842497 > 2^59. | |
409 | ||
410 | The overall period exceeds 2^123. | |
411 | ||
412 | http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu | |
413 | ||
414 | The above web site has an archive of a newsgroup posting from George | |
415 | Marsaglia with the statement: | |
416 | ||
417 | Subject: Random numbers for C: Improvements. | |
418 | Date: Fri, 15 Jan 1999 11:41:47 -0500 | |
419 | From: George Marsaglia <geo@stat.fsu.edu> | |
420 | Message-ID: <369F6FCA.74C7C041@stat.fsu.edu> | |
421 | References: <369B5E30.65A55FD1@stat.fsu.edu> | |
422 | Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis | |
423 | Lines: 93 | |
424 | ||
425 | As I hoped, several suggestions have led to | |
426 | improvements in the code for RNG's I proposed for | |
427 | use in C. (See the thread "Random numbers for C: Some | |
428 | suggestions" in previous postings.) The improved code | |
429 | is listed below. | |
430 | ||
431 | A question of copyright has also been raised. Unlike | |
432 | DIEHARD, there is no copyright on the code below. You | |
433 | are free to use it in any way you want, but you may | |
434 | wish to acknowledge the source, as a courtesy. | |
435 | ||
436 | "There is no copyright on the code below." included the original | |
cdaa9fc4 | 437 | KISS algorithm. */ |
5f251c26 | 438 | |
5f251c26 SK |
439 | #define GFC_SL(k, n) ((k)^((k)<<(n))) |
440 | #define GFC_SR(k, n) ((k)^((k)>>(n))) | |
441 | ||
442 | static const GFC_INTEGER_4 kiss_size = 4; | |
4b6903ec | 443 | #define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069} |
5f251c26 SK |
444 | static const GFC_UINTEGER_4 kiss_default_seed[4] = KISS_DEFAULT_SEED; |
445 | static GFC_UINTEGER_4 kiss_seed[4] = KISS_DEFAULT_SEED; | |
446 | ||
447 | /* kiss_random_kernel() returns an integer value in the range of | |
448 | (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers | |
449 | should be uniform. */ | |
450 | ||
451 | static GFC_UINTEGER_4 | |
452 | kiss_random_kernel(void) | |
453 | { | |
5f251c26 SK |
454 | GFC_UINTEGER_4 kiss; |
455 | ||
456 | kiss_seed[0] = 69069 * kiss_seed[0] + 1327217885; | |
457 | kiss_seed[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed[1],13),17),5); | |
458 | kiss_seed[2] = 18000 * (kiss_seed[2] & 65535) + (kiss_seed[2] >> 16); | |
459 | kiss_seed[3] = 30903 * (kiss_seed[3] & 65535) + (kiss_seed[3] >> 16); | |
460 | kiss = kiss_seed[0] + kiss_seed[1] + (kiss_seed[2] << 16) + kiss_seed[3]; | |
461 | ||
462 | return kiss; | |
5f251c26 SK |
463 | } |
464 | ||
cdaa9fc4 | 465 | /* This function produces a REAL(4) value from the uniform distribution |
5f251c26 SK |
466 | with range [0,1). */ |
467 | ||
468 | void | |
7d7b8bfe | 469 | random_r4 (GFC_REAL_4 *x) |
5f251c26 | 470 | { |
5f251c26 SK |
471 | GFC_UINTEGER_4 kiss; |
472 | ||
a9e7b9d3 PB |
473 | kiss = kiss_random_kernel (); |
474 | /* Burn a random number, so the REAL*4 and REAL*8 functions | |
475 | produce similar sequences of random numbers. */ | |
476 | kiss_random_kernel (); | |
477 | *x = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0); | |
5f251c26 | 478 | } |
7d7b8bfe | 479 | iexport(random_r4); |
5f251c26 SK |
480 | |
481 | /* This function produces a REAL(8) value from the uniform distribution | |
482 | with range [0,1). */ | |
483 | ||
484 | void | |
7d7b8bfe | 485 | random_r8 (GFC_REAL_8 *x) |
5f251c26 | 486 | { |
5f251c26 SK |
487 | GFC_UINTEGER_8 kiss; |
488 | ||
a9e7b9d3 PB |
489 | kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32; |
490 | kiss += kiss_random_kernel (); | |
491 | *x = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0); | |
5f251c26 | 492 | } |
7d7b8bfe | 493 | iexport(random_r8); |
5f251c26 SK |
494 | |
495 | /* This function fills a REAL(4) array with values from the uniform | |
496 | distribution with range [0,1). */ | |
497 | ||
498 | void | |
7d7b8bfe | 499 | arandom_r4 (gfc_array_r4 *x) |
5f251c26 | 500 | { |
e33e218b TK |
501 | index_type count[GFC_MAX_DIMENSIONS]; |
502 | index_type extent[GFC_MAX_DIMENSIONS]; | |
503 | index_type stride[GFC_MAX_DIMENSIONS]; | |
5f251c26 SK |
504 | index_type stride0; |
505 | index_type dim; | |
506 | GFC_REAL_4 *dest; | |
507 | int n; | |
508 | ||
509 | dest = x->data; | |
510 | ||
511 | if (x->dim[0].stride == 0) | |
512 | x->dim[0].stride = 1; | |
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 | while (dest) | |
528 | { | |
7d7b8bfe | 529 | random_r4 (dest); |
5f251c26 SK |
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 | } | |
558 | ||
cdaa9fc4 | 559 | /* This function fills a REAL(8) array with values from the uniform |
5f251c26 SK |
560 | distribution with range [0,1). */ |
561 | ||
562 | void | |
7d7b8bfe | 563 | arandom_r8 (gfc_array_r8 *x) |
5f251c26 | 564 | { |
e33e218b TK |
565 | index_type count[GFC_MAX_DIMENSIONS]; |
566 | index_type extent[GFC_MAX_DIMENSIONS]; | |
567 | index_type stride[GFC_MAX_DIMENSIONS]; | |
5f251c26 SK |
568 | index_type stride0; |
569 | index_type dim; | |
570 | GFC_REAL_8 *dest; | |
571 | int n; | |
572 | ||
573 | dest = x->data; | |
574 | ||
575 | if (x->dim[0].stride == 0) | |
576 | x->dim[0].stride = 1; | |
577 | ||
578 | dim = GFC_DESCRIPTOR_RANK (x); | |
579 | ||
580 | for (n = 0; n < dim; n++) | |
581 | { | |
582 | count[n] = 0; | |
583 | stride[n] = x->dim[n].stride; | |
584 | extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound; | |
585 | if (extent[n] <= 0) | |
586 | return; | |
587 | } | |
588 | ||
589 | stride0 = stride[0]; | |
590 | ||
591 | while (dest) | |
592 | { | |
7d7b8bfe | 593 | random_r8 (dest); |
5f251c26 SK |
594 | |
595 | /* Advance to the next element. */ | |
596 | dest += stride0; | |
597 | count[0]++; | |
598 | /* Advance to the next source element. */ | |
599 | n = 0; | |
600 | while (count[n] == extent[n]) | |
601 | { | |
602 | /* When we get to the end of a dimension, reset it and increment | |
603 | the next dimension. */ | |
604 | count[n] = 0; | |
605 | /* We could precalculate these products, but this is a less | |
606 | frequently used path so probably not worth it. */ | |
607 | dest -= stride[n] * extent[n]; | |
608 | n++; | |
609 | if (n == dim) | |
610 | { | |
611 | dest = NULL; | |
612 | break; | |
613 | } | |
614 | else | |
615 | { | |
616 | count[n]++; | |
617 | dest += stride[n]; | |
618 | } | |
619 | } | |
620 | } | |
621 | } | |
622 | ||
7d7b8bfe | 623 | /* random_seed is used to seed the PRNG with either a default |
420aa7b8 | 624 | set of seeds or user specified set of seeds. random_seed |
5f251c26 SK |
625 | must be called with no argument or exactly one argument. */ |
626 | ||
627 | void | |
f21edfd6 | 628 | random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get) |
5f251c26 | 629 | { |
5f251c26 SK |
630 | int i; |
631 | ||
632 | if (size == NULL && put == NULL && get == NULL) | |
633 | { | |
634 | /* From the standard: "If no argument is present, the processor assigns | |
cdaa9fc4 | 635 | a processor-dependent value to the seed." */ |
5f251c26 SK |
636 | kiss_seed[0] = kiss_default_seed[0]; |
637 | kiss_seed[1] = kiss_default_seed[1]; | |
638 | kiss_seed[2] = kiss_default_seed[2]; | |
639 | kiss_seed[3] = kiss_default_seed[3]; | |
640 | } | |
641 | ||
642 | if (size != NULL) | |
643 | *size = kiss_size; | |
644 | ||
645 | if (put != NULL) | |
646 | { | |
cdaa9fc4 | 647 | /* If the rank of the array is not 1, abort. */ |
5f251c26 SK |
648 | if (GFC_DESCRIPTOR_RANK (put) != 1) |
649 | runtime_error ("Array rank of PUT is not 1."); | |
650 | ||
cdaa9fc4 | 651 | /* If the array is too small, abort. */ |
5f251c26 SK |
652 | if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size) |
653 | runtime_error ("Array size of PUT is too small."); | |
654 | ||
655 | if (put->dim[0].stride == 0) | |
656 | put->dim[0].stride = 1; | |
657 | ||
cdaa9fc4 | 658 | /* This code now should do correct strides. */ |
5f251c26 SK |
659 | for (i = 0; i < kiss_size; i++) |
660 | kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride]; | |
661 | } | |
662 | ||
cdaa9fc4 | 663 | /* Return the seed to GET data. */ |
5f251c26 SK |
664 | if (get != NULL) |
665 | { | |
cdaa9fc4 | 666 | /* If the rank of the array is not 1, abort. */ |
5f251c26 SK |
667 | if (GFC_DESCRIPTOR_RANK (get) != 1) |
668 | runtime_error ("Array rank of GET is not 1."); | |
669 | ||
cdaa9fc4 | 670 | /* If the array is too small, abort. */ |
5f251c26 SK |
671 | if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size) |
672 | runtime_error ("Array size of GET is too small."); | |
673 | ||
674 | if (get->dim[0].stride == 0) | |
675 | get->dim[0].stride = 1; | |
676 | ||
cdaa9fc4 | 677 | /* This code now should do correct strides. */ |
5f251c26 SK |
678 | for (i = 0; i < kiss_size; i++) |
679 | get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i]; | |
680 | } | |
681 | } | |
7d7b8bfe RH |
682 | iexport(random_seed); |
683 | ||
684 | #endif /* mersenne twister */ |