OSDN Git Service

fdb99cc040ef092046d10904a4bd0219da20dfe6
[pf3gnuchains/gcc-fork.git] / libgfortran / intrinsics / random.c
1 /* Implementation of the RANDOM intrinsics
2    Copyright 2002, 2004 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., 59 Temple Place - Suite 330,
30 Boston, MA 02111-1307, USA.  */
31
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
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
71 #include <stdio.h>
72 #include <stdlib.h>
73 #include <sys/types.h>
74 #include <sys/stat.h>
75 #include <fcntl.h>
76
77 #ifdef HAVE_UNISTD_H
78 #include <unistd.h>
79 #endif
80
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
111 void
112 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
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 */
161       if (put->dim[0].stride == 0)
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 }
189 iexport(random_seed);
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
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 }
231 iexport(random_r4);
232
233 /* A routine to return a REAL(KIND=8).  */
234
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 }
247 iexport(random_r8);
248
249 /* Code to handle arrays will follow here.  */
250
251 /* REAL(KIND=4) REAL array.  */
252
253 void
254 arandom_r4 (gfc_array_r4 * harv)
255 {
256   index_type count[GFC_MAX_DIMENSIONS - 1];
257   index_type extent[GFC_MAX_DIMENSIONS - 1];
258   index_type stride[GFC_MAX_DIMENSIONS - 1];
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
326 void
327 arandom_r8 (gfc_array_r8 * harv)
328 {
329   index_type count[GFC_MAX_DIMENSIONS - 1];
330   index_type extent[GFC_MAX_DIMENSIONS - 1];
331   index_type stride[GFC_MAX_DIMENSIONS - 1];
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 }
397
398 #else
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
437 KISS algorithm.  */
438
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;
443 #define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069};
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 {
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;
463 }
464
465 /*  This function produces a REAL(4) value from the uniform distribution
466     with range [0,1).  */
467
468 void
469 random_r4 (GFC_REAL_4 *x)
470 {
471   GFC_UINTEGER_4 kiss;
472
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);
478 }
479 iexport(random_r4);
480
481 /*  This function produces a REAL(8) value from the uniform distribution
482     with range [0,1).  */
483
484 void
485 random_r8 (GFC_REAL_8 *x)
486 {
487   GFC_UINTEGER_8 kiss;
488
489   kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
490   kiss += kiss_random_kernel ();
491   *x = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
492 }
493 iexport(random_r8);
494
495 /*  This function fills a REAL(4) array with values from the uniform
496     distribution with range [0,1).  */
497
498 void
499 arandom_r4 (gfc_array_r4 *x)
500 {
501   index_type count[GFC_MAX_DIMENSIONS - 1];
502   index_type extent[GFC_MAX_DIMENSIONS - 1];
503   index_type stride[GFC_MAX_DIMENSIONS - 1];
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     {
529       random_r4 (dest);
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
559 /*  This function fills a REAL(8) array with values from the uniform
560     distribution with range [0,1).  */
561
562 void
563 arandom_r8 (gfc_array_r8 *x)
564 {
565   index_type count[GFC_MAX_DIMENSIONS - 1];
566   index_type extent[GFC_MAX_DIMENSIONS - 1];
567   index_type stride[GFC_MAX_DIMENSIONS - 1];
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     {
593       random_r8 (dest);
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
623 /* random_seed is used to seed the PRNG with either a default
624    set of seeds or user specified set of seeds.  random_seed 
625    must be called with no argument or exactly one argument.  */
626
627 void
628 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
629 {
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
635          a processor-dependent value to the seed."  */
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     {
647       /* If the rank of the array is not 1, abort.  */
648       if (GFC_DESCRIPTOR_RANK (put) != 1)
649         runtime_error ("Array rank of PUT is not 1.");
650
651       /* If the array is too small, abort.  */
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
658       /*  This code now should do correct strides.  */
659       for (i = 0; i < kiss_size; i++)
660         kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
661     }
662
663   /* Return the seed to GET data.  */
664   if (get != NULL)
665     {
666       /* If the rank of the array is not 1, abort.  */
667       if (GFC_DESCRIPTOR_RANK (get) != 1)
668         runtime_error ("Array rank of GET is not 1.");
669
670       /* If the array is too small, abort.  */
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
677       /*  This code now should do correct strides.  */
678       for (i = 0; i < kiss_size; i++)
679         get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
680     }
681 }
682 iexport(random_seed);
683
684 #endif /* mersenne twister */