This is the mail archive of the gcc-patches@gcc.gnu.org mailing list for the GCC project.
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
| Other format: | [Raw text] | |
> See attached.
> fortran.diff is applied to gcc/fortran.
> random.c.diff is applied to random.c
I tweaked the int->float conversion routines, and applied as follows.
Paul
2004-05-30 Steven G. Kargl <kargls@comcast.net>
* iresolve.c (gfc_resolve_random_number): Clean up conditional.
libgfortran/
* libgfortran.h (random_seed): Update prototype.
* intrinsics/random.c: Disable old implementation and add new one.
testsuite/
* gfortran.fortran-torture/execute/random_1.f90: New test.
Index: gcc/fortran/iresolve.c
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/gcc/fortran/iresolve.c,v
retrieving revision 1.5
diff -u -p -r1.5 iresolve.c
--- gcc/fortran/iresolve.c 22 May 2004 12:47:38 -0000 1.5
+++ gcc/fortran/iresolve.c 30 May 2004 10:14:04 -0000
@@ -1363,12 +1363,15 @@ gfc_resolve_random_number (gfc_code * c
int kind;
kind = c->ext.actual->expr->ts.kind;
- name = gfc_get_string ((c->ext.actual->expr->rank == 0) ?
- PREFIX("random_r%d") : PREFIX("arandom_r%d"),
- kind);
+ if (c->ext.actual->expr->rank == 0)
+ name = gfc_get_string (PREFIX("random_r%d"), kind);
+ else
+ name = gfc_get_string (PREFIX("arandom_r%d"), kind);
+
c->resolved_sym = gfc_get_intrinsic_sub_symbol (name);
}
+
/* Determine if the arguments to SYSTEM_CLOCK are INTEGER(4) or INTEGER(8) */
void
Index: libgfortran/libgfortran.h
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/libgfortran.h,v
retrieving revision 1.4
diff -u -p -r1.4 libgfortran.h
--- libgfortran/libgfortran.h 18 May 2004 16:06:08 -0000 1.4
+++ libgfortran/libgfortran.h 30 May 2004 10:32:26 -0000
@@ -399,8 +399,8 @@ GFC_INTEGER_4 compare_string (GFC_INTEGE
/* random.c */
#define random_seed prefix(random_seed)
-void random_seed (GFC_INTEGER_4 * size, const gfc_array_i4 * put,
- const gfc_array_i4 * get);
+void random_seed (GFC_INTEGER_4 * size, gfc_array_i4 * put,
+ gfc_array_i4 * get);
#endif
Index: libgfortran/intrinsics/random.c
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/intrinsics/random.c,v
retrieving revision 1.5
diff -u -p -r1.5 random.c
--- libgfortran/intrinsics/random.c 30 May 2004 09:53:09 -0000 1.5
+++ libgfortran/intrinsics/random.c 30 May 2004 10:42:50 -0000
@@ -1,18 +1,7 @@
/* Implementation of the RANDOM intrinsics
Copyright 2002, 2004 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 ).
+ and Steve Kargl.
This file is part of the GNU Fortran 95 runtime library (libgfortran).
@@ -31,6 +20,31 @@ License along with libgfor; see the file
write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA. */
+#if 0
+
+/* The Mersenne Twister code is currently commented out due to
+
+ (1) Simple user specified seeds lead to really bad sequences for
+ nearly 100000 random numbers.
+ (2) open(), read(), and close() are not properly declared via header
+ files.
+ (3) The global index i is abused and causes unexpected behavior with
+ GET and PUT.
+ (4) See PR 15619.
+
+ 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>
@@ -362,4 +376,306 @@ arandom_r8 (gfc_array_r8 * harv)
}
}
}
+#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.
+
+"There is no copyright on the code below." included the original
+KISS algorithm. */
+
+#include "config.h"
+#include "libgfortran.h"
+
+#define GFC_SL(k, n) ((k)^((k)<<(n)))
+#define GFC_SR(k, n) ((k)^((k)>>(n)))
+
+static const GFC_INTEGER_4 kiss_size = 4;
+#define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069};
+static const GFC_UINTEGER_4 kiss_default_seed[4] = KISS_DEFAULT_SEED;
+static GFC_UINTEGER_4 kiss_seed[4] = KISS_DEFAULT_SEED;
+
+/* kiss_random_kernel() returns an integer value in the range of
+ (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
+ should be uniform. */
+
+static GFC_UINTEGER_4
+kiss_random_kernel(void)
+{
+
+ GFC_UINTEGER_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)
+{
+
+ GFC_UINTEGER_4 kiss;
+
+ do
+ {
+ kiss = kiss_random_kernel ();
+ *x = (GFC_REAL_4)kiss / (GFC_REAL_4)(~(GFC_UINTEGER_4) 0);
+ }
+ while (*x == 1.0);
+
+}
+
+/* This function produces a REAL(8) value from the uniform distribution
+ with range [0,1). */
+
+void
+prefix(random_r8) (GFC_REAL_8 *x)
+{
+
+ GFC_UINTEGER_8 kiss;
+
+ do
+ {
+ kiss = (((GFC_UINTEGER_8)kiss_random_kernel ()) << 32)
+ + kiss_random_kernel ();
+ *x = (GFC_REAL_8)kiss / (GFC_REAL_8)(~(GFC_UINTEGER_8) 0);
+ }
+ while (*x != 0);
+
+}
+
+/* This function fills a REAL(4) array with values 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
+random_seed (GFC_INTEGER_4 *size, gfc_array_i4 * put,
+ 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)
+ runtime_error ("Array rank of PUT is not 1.");
+
+ /* If the array is too small, abort */
+ if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
+ runtime_error ("Array size of PUT is too small.");
+
+ if (put->dim[0].stride == 0)
+ put->dim[0].stride = 1;
+
+ /* This code now should do correct strides. */
+ for (i = 0; i < kiss_size; i++)
+ kiss_seed[i] =(GFC_UINTEGER_4) 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)
+ runtime_error ("Array rank of GET is not 1.");
+
+ /* If the array is too small, abort. */
+ if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
+ runtime_error ("Array size of GET is too small.");
+
+ if (get->dim[0].stride == 0)
+ get->dim[0].stride = 1;
+
+ /* This code now should do correct strides. */
+ for (i = 0; i < kiss_size; i++)
+ get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
+ }
+}
+
Attachment:
random_1.f90
Description: Text document
| Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
|---|---|---|
| Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |