[gfortran] PATCH -- Fix for PR 15619
Steve Kargl
sgk@troutmask.apl.washington.edu
Thu May 27 10:05:00 GMT 2004
On Thu, May 27, 2004 at 12:12:59AM +0100, Paul Brook wrote:
>
> No, still wrong. Fails the test appended to this mail.
>
> Possible solution: add an extra word to the seed which contains the current
> position. Obviously you'll need to cope with arbitary values when the user
> specifies PUT=, but curent_index=abs(current_index) % N should do the trick.
>
> Oh, and while you're at it, please rename the global variable "i" to
> something more sane ;)
>
Sigh. I really didn't want to become an expert on PRNG. I only
want to use them. :-) So, here's my current solution.
See the attached diff. It replaces the Mersenne Twister with
Marsaglia's KISS PRNG. I left Lars' code in random.c where I
simply enclosed it in an #if 0 ... #endif. KISS seems to fix
PR 15619, the code I elided from this email, and the seeding
problem I discuss in other emails. I've retained Lars'
function names of prefix(random_r4), prefix(arandom_r4),
prefix(random_r8), and prefix(arandom_r8); although I think
the names should be more descriptivei, e.g., prefix(random_number_r4)
and prefix(random_number_array_r4). This will allow Lars or somone
else to continue hacking on the Mersenne Twister without changing the API.
2004-05-26 Steven G. Kargl <kargls@comcast.net>
* random.c: Replace Mersenne Twister with Marsaglia's KISS PRNG
(kiss_random_kernel): New function.
(random_r4,random_r8): Use it.
(random_seed,arandom_r4,arandom_r8): Update to KISS
--
Steve
-------------- next part --------------
Index: random.c
===================================================================
RCS file: /cvs/gcc/gcc/libgfortran/intrinsics/random.c,v
retrieving revision 1.4
diff -u -r1.4 random.c
--- random.c 23 May 2004 16:18:22 -0000 1.4
+++ random.c 27 May 2004 00:19:33 -0000
@@ -2,17 +2,6 @@
Copyright 2002 Free Software Foundation, Inc.
Contributed by Lars Segerlund <seger@linuxmail.org>
- The algorithm was taken from the paper :
-
- Mersenne Twister: 623-dimensionally equidistributed
- uniform pseudorandom generator.
-
- by: Makoto Matsumoto
- Takuji Nishimura
-
- Which appeared in the: ACM Transactions on Modelling and Computer
- Simulations: Special Issue on Uniform Random Number
- Generation. ( Early in 1998 ).
This file is part of the GNU Fortran 95 runtime library (libgfortran).
@@ -31,10 +20,27 @@
write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA. */
+#if 0
+/* The algorithm was taken from the paper :
+
+ Mersenne Twister: 623-dimensionally equidistributed
+ uniform pseudorandom generator.
+
+ by: Makoto Matsumoto
+ Takuji Nishimura
+
+ Which appeared in the: ACM Transactions on Modelling and Computer
+ Simulations: Special Issue on Uniform Random Number
+ Generation. ( Early in 1998 ). */
+
#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
+#ifdef __FreeBSD__
+# include <sys/uio.h>
+# include <unistd.h>
+#endif
#include <sys/stat.h>
#include <fcntl.h>
#include <assert.h>
@@ -74,6 +80,9 @@
random_seed (GFC_INTEGER_4 * size, const gfc_array_i4 * put,
const gfc_array_i4 * get)
{
+
+ int j;
+
/* Initialize the seed in system dependent manner. */
if (get == NULL && put == NULL && size == NULL)
{
@@ -83,10 +92,10 @@
{
/* We dont have urandom. */
GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed;
- for (i = 0; i < N; i++)
+ for (j = 0; j < N; j++)
{
s = s * 29943829 - 1;
- seed[i] = s;
+ seed[j] = s;
}
}
else
@@ -125,8 +134,10 @@
return;
/* This code now should do correct strides. */
- for (i = 0; i < N; i++)
- seed[i] = put->data[i * put->dim[0].stride];
+ for (j = 0; j < N; j++)
+ seed[j] = put->data[j * put->dim[0].stride];
+
+ i = N;
}
/* Return the seed to GET data */
@@ -145,8 +156,8 @@
return;
/* This code now should do correct strides. */
- for (i = 0; i < N; i++)
- get->data[i * get->dim[0].stride] = seed[i];
+ for (j = 0; j < N; j++)
+ get->data[j * get->dim[0].stride] = seed[j];
}
}
@@ -186,6 +197,7 @@
/* Regenerate if we need to. */
if (i >= N)
random_generate ();
+fprintf(stderr, "i = %d\n", i);
/* Convert uint32 to REAL(KIND=4). */
*harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
@@ -359,3 +371,287 @@
}
}
+#endif /* Mersenne Twister code */
+
+/* George Marsaglia's KISS (Keep It Simple Stupid) random number generator.
+
+ This PRNG combines:
+
+ (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
+ of 2^32,
+ (2) A 3-shift shift-register generator with a period of 2^32-1,
+ (3) Two 16-bit multiply-with-carry generators with a period of
+ 597273182964842497 > 2^59.
+
+ The overall period exceeds 2^123.
+
+ http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
+
+ The above web site has an archive of a newsgroup posting from George
+ Marsaglia with the statement:
+
+ Subject: Random numbers for C: Improvements.
+ Date: Fri, 15 Jan 1999 11:41:47 -0500
+ From: George Marsaglia <geo@stat.fsu.edu>
+ Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
+ References: <369B5E30.65A55FD1@stat.fsu.edu>
+ Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
+ Lines: 93
+
+ As I hoped, several suggestions have led to
+ improvements in the code for RNG's I proposed for
+ use in C. (See the thread "Random numbers for C: Some
+ suggestions" in previous postings.) The improved code
+ is listed below.
+
+ A question of copyright has also been raised. Unlike
+ DIEHARD, there is no copyright on the code below. You
+ are free to use it in any way you want, but you may
+ wish to acknowledge the source, as a courtesy.
+*/
+#include "config.h"
+#include <stdlib.h> /*Needed for abort() */
+#include "libgfortran.h"
+
+#define GFC_SL(k, n) ((k)^((k)<<(n)))
+#define GFC_SR(k, n) ((k)^((k)>>(n)))
+
+const GFC_INTEGER_4 kiss_size = 4;
+const GFC_INTEGER_4
+ kiss_default_seed[4] = {123456789, 362436069, 521288629, 916191069};
+static GFC_INTEGER_4
+ kiss_seed[4] = {123456789, 362436069, 521288629, 916191069};
+
+/* kiss_random_kernel() returns an integer value in the range of
+ [- GFC_INTEGER_4_HUGE, GFC_INTEGER_4_HUGE]. The distribution of
+ pseudorandom numbers should be uniform. */
+
+static GFC_INTEGER_4 kiss_random_kernel(void)
+{
+
+ GFC_INTEGER_4 kiss;
+
+ kiss_seed[0] = 69069 * kiss_seed[0] + 1327217885;
+ kiss_seed[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed[1],13),17),5);
+ kiss_seed[2] = 18000 * (kiss_seed[2] & 65535) + (kiss_seed[2] >> 16);
+ kiss_seed[3] = 30903 * (kiss_seed[3] & 65535) + (kiss_seed[3] >> 16);
+ kiss = kiss_seed[0] + kiss_seed[1] + (kiss_seed[2] << 16) + kiss_seed[3];
+
+ return kiss;
+
+}
+
+/* This function produces a REAL(4) value in the uniform distribution
+ with range [0,1). */
+
+void prefix(random_r4) (GFC_REAL_4 *x)
+{
+
+ *x = (GFC_REAL_4) kiss_random_kernel();
+ *x = ((*x < 0) ? - *x : *x) / ((GFC_REAL_4) (GFC_INTEGER_4_HUGE) + 1.);
+
+}
+
+/* This function produces a REAL(8) value from the uniform distribution
+ with range [0,1). */
+
+void prefix(random_r8) (GFC_REAL_8 *x)
+{
+
+ *x = (GFC_REAL_8) kiss_random_kernel();
+ *x = ((*x < 0) ? - *x : *x) / ((GFC_REAL_8) (GFC_INTEGER_4_HUGE) + 1.);
+
+}
+
+/* This function fills a REAL(4) array with valuse from the uniform
+ distribution with range [0,1). */
+
+void prefix(arandom_r4) (gfc_array_r4 *x)
+{
+
+ index_type count[GFC_MAX_DIMENSIONS - 1];
+ index_type extent[GFC_MAX_DIMENSIONS - 1];
+ index_type stride[GFC_MAX_DIMENSIONS - 1];
+ index_type stride0;
+ index_type dim;
+ GFC_REAL_4 *dest;
+ int n;
+
+ dest = x->data;
+
+ if (x->dim[0].stride == 0)
+ x->dim[0].stride = 1;
+
+ dim = GFC_DESCRIPTOR_RANK (x);
+
+ for (n = 0; n < dim; n++)
+ {
+ count[n] = 0;
+ stride[n] = x->dim[n].stride;
+ extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
+ if (extent[n] <= 0)
+ return;
+ }
+
+ stride0 = stride[0];
+
+ while (dest)
+ {
+ prefix(random_r4) (dest);
+
+ /* Advance to the next element. */
+ dest += stride0;
+ count[0]++;
+ /* Advance to the next source element. */
+ n = 0;
+ while (count[n] == extent[n])
+ {
+ /* When we get to the end of a dimension, reset it and increment
+ the next dimension. */
+ count[n] = 0;
+ /* We could precalculate these products, but this is a less
+ frequently used path so probably not worth it. */
+ dest -= stride[n] * extent[n];
+ n++;
+ if (n == dim)
+ {
+ dest = NULL;
+ break;
+ }
+ else
+ {
+ count[n]++;
+ dest += stride[n];
+ }
+ }
+ }
+}
+
+/* This function fills a REAL(8) array with valuse from the uniform
+ distribution with range [0,1). */
+
+void prefix(arandom_r8) (gfc_array_r8 *x)
+{
+
+ index_type count[GFC_MAX_DIMENSIONS - 1];
+ index_type extent[GFC_MAX_DIMENSIONS - 1];
+ index_type stride[GFC_MAX_DIMENSIONS - 1];
+ index_type stride0;
+ index_type dim;
+ GFC_REAL_8 *dest;
+ int n;
+
+ dest = x->data;
+
+ if (x->dim[0].stride == 0)
+ x->dim[0].stride = 1;
+
+ dim = GFC_DESCRIPTOR_RANK (x);
+
+ for (n = 0; n < dim; n++)
+ {
+ count[n] = 0;
+ stride[n] = x->dim[n].stride;
+ extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
+ if (extent[n] <= 0)
+ return;
+ }
+
+ stride0 = stride[0];
+
+ while (dest)
+ {
+ prefix(random_r8) (dest);
+
+ /* Advance to the next element. */
+ dest += stride0;
+ count[0]++;
+ /* Advance to the next source element. */
+ n = 0;
+ while (count[n] == extent[n])
+ {
+ /* When we get to the end of a dimension, reset it and increment
+ the next dimension. */
+ count[n] = 0;
+ /* We could precalculate these products, but this is a less
+ frequently used path so probably not worth it. */
+ dest -= stride[n] * extent[n];
+ n++;
+ if (n == dim)
+ {
+ dest = NULL;
+ break;
+ }
+ else
+ {
+ count[n]++;
+ dest += stride[n];
+ }
+ }
+ }
+}
+
+/* prefix(random_seed) is used to seed the PRNG with either a default
+ set of seeds or user specified set of seed. prefix(random_seed)
+ must be called with no argument or exactly one argument. */
+
+void prefix(random_seed) (GFC_INTEGER_4 *size, const gfc_array_i4 * put,
+ const gfc_array_i4 * get)
+{
+
+ int i;
+
+ if (size == NULL && put == NULL && get == NULL)
+ {
+ /* From the standard: "If no argument is present, the processor assigns
+ a processor-dependent value to the seed." */
+ kiss_seed[0] = kiss_default_seed[0];
+ kiss_seed[1] = kiss_default_seed[1];
+ kiss_seed[2] = kiss_default_seed[2];
+ kiss_seed[3] = kiss_default_seed[3];
+ }
+
+ if (size != NULL)
+ *size = kiss_size;
+
+ if (put != NULL)
+ {
+ /* If the rank of the array is not 1, abort */
+ if (GFC_DESCRIPTOR_RANK (put) != 1)
+ abort ();
+
+ /* If the array is too small, abort */
+ if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
+ abort ();
+
+ /* If this is the case, the array is a temporary.
+ XXX Is this correct? */
+ if (put->dim[0].stride == 0)
+ return;
+
+ /* This code now should do correct strides. */
+ for (i = 0; i < kiss_size; i++)
+ kiss_seed[i] = put->data[i * put->dim[0].stride];
+ }
+
+ /* Return the seed to GET data */
+ if (get != NULL)
+ {
+ /* If the rank of the array is not 1, abort. */
+ if (GFC_DESCRIPTOR_RANK (get) != 1)
+ abort ();
+
+ /* If the array is too small, abort. */
+ if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
+ abort ();
+
+ /* If this is the case, the array is a temporary.
+ XXX Is this correct? */
+ if (get->dim[0].stride == 0)
+ return;
+
+ /* This code now should do correct strides. */
+ for (i = 0; i < kiss_size; i++)
+ get->data[i * get->dim[0].stride] = kiss_seed[i];
+ }
+}
More information about the Gcc-patches
mailing list