Re: [gfortran] PATCH -- random.c: fix infinite loop and make sequences consistent

[Paul Brook <paul at codesourcery dot com>]

On Sunday 13 June 2004 18:16, Steven G. Kargl wrote:
> Gang,
>
> The attached patch fixes an infinite loop in the random_number
> routine for REAL*8 arrays and it burns a random number in
> the REAL*4 routine to keep the REAL*4 and REAL*8 sequence consistent.
>
>
> 2004-06-13  Steven G. Kargl  <kargls@comcast.net>
>      * random.c (random_r4): Burn a random number
>        (random_r8): fix infinite loop.

Still not quite there I'm aftraid. It's possible to generate vales where 
(float)x ==1.0 && (double)x<1.0. This will cause the sequences to diverge.

I've implemented a int->float conversion routine along the lines suggested by 
Richard Henderson.

I also added a testcase.

Tested on i686-linux. Applied to mainline.

Paul

2004-06-13  Paul Brook  <paul@codesourcery.com>

	* Makefile.am (gfor_helper_src): Add runtime/normalize.f90.
	* configure.ac: Add checks for nextafter and nextafterf.
	* Makefile.in, config.h.in, configure: Regenerate.
	* libgfortran.h (normalize_r4_i4, normalize_r8_i8): Declare.
	* intrinsics/rand.c (rand): Use normalize_r4_i4.
	* intrinsics/random.c (random_r4): Use normalize_r4_i4.
	(random_r8): Use normalize_r8_i8.
	* runtime/normalize.c: New file.
testsuite/
	* gfortran.fortran-torture/execute/random_2.f90: New test.

Index: gcc/testsuite/gfortran.fortran-torture/execute/random_2.f90
===================================================================
RCS file: gcc/testsuite/gfortran.fortran-torture/execute/random_2.f90
diff -N gcc/testsuite/gfortran.fortran-torture/execute/random_2.f90
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ gcc/testsuite/gfortran.fortran-torture/execute/random_2.f90	13 Jun 2004 
20:51:32 -0000
@@ -0,0 +1,24 @@
+! Check that the real(4) and real(8) random number generators return the same
+! sequence of values.
+program random_4
+  integer, dimension(:), allocatable :: seed
+  real(kind=4), dimension(10) :: r4
+  real(kind=8), dimension(10) :: r8
+  real, parameter :: delta = 0.0001
+  integer n
+
+  call random_seed (size=n)
+  allocate (seed(n))
+  call random_seed (get=seed)
+  ! Test both array valued and scalar routines.
+  call random_number(r4)
+  call random_number (r4(10))
+
+  ! Reset the seed and get the real(8) values.
+  call random_seed (put=seed)
+  call random_number(r8)
+  call random_number (r8(10))
+
+  if (any ((r4 - r8) .gt. delta)) call abort
+end program
+
Index: libgfortran/Makefile.am
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/Makefile.am,v
retrieving revision 1.10
diff -u -p -r1.10 Makefile.am
--- libgfortran/Makefile.am	12 Jun 2004 17:59:29 -0000	1.10
+++ libgfortran/Makefile.am	13 Jun 2004 22:23:34 -0000
@@ -58,7 +58,8 @@ intrinsics/system_clock.c \
 intrinsics/transpose_generic.c \
 intrinsics/unpack_generic.c \
 runtime/in_pack_generic.c \
-runtime/in_unpack_generic.c
+runtime/in_unpack_generic.c \
+runtime/normalize.c
 
 gfor_src= \
 runtime/environ.c \
Index: libgfortran/configure.ac
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/configure.ac,v
retrieving revision 1.2
diff -u -p -r1.2 configure.ac
--- libgfortran/configure.ac	12 Jun 2004 17:59:30 -0000	1.2
+++ libgfortran/configure.ac	13 Jun 2004 20:08:31 -0000
@@ -169,6 +169,9 @@ AC_CHECK_FUNCS(getrusage times)
 # Check for some C99 functions
 AC_CHECK_LIB([m],[round],[AC_DEFINE([HAVE_ROUND],[1],["c99 function"])])
 AC_CHECK_LIB([m],[roundf],[AC_DEFINE([HAVE_ROUNDF],[1],["c99 function"])])
+# And other IEEE math functions
+AC_CHECK_LIB([m],[nextafter],[AC_DEFINE([HAVE_NEXTAFTER],[1],[libm includes 
nextafter])])
+AC_CHECK_LIB([m],[nextafterf],[AC_DEFINE([HAVE_NEXTAFTERF],[1],[libm includes 
nextafterf])])
 
 # Let the user override this
 AC_ARG_ENABLE(cmath,
Index: libgfortran/libgfortran.h
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/libgfortran.h,v
retrieving revision 1.6
diff -u -p -r1.6 libgfortran.h
--- libgfortran/libgfortran.h	12 Jun 2004 15:15:41 -0000	1.6
+++ libgfortran/libgfortran.h	13 Jun 2004 20:05:47 -0000
@@ -408,5 +408,13 @@ GFC_INTEGER_4 compare_string (GFC_INTEGE
 void random_seed (GFC_INTEGER_4 * size, gfc_array_i4 * put,
 		  gfc_array_i4 * get);
 
+/* normalize.c */
+
+#define normalize_r4_i4 prefix(normalize_r4_i4)
+GFC_REAL_4 normalize_r4_i4 (GFC_UINTEGER_4, GFC_UINTEGER_4);
+
+#define normalize_r8_i8 prefix(normalize_r8_i8)
+GFC_REAL_8 normalize_r8_i8 (GFC_UINTEGER_8, GFC_UINTEGER_8);
+
 #endif
 
Index: libgfortran/intrinsics/rand.c
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/intrinsics/rand.c,v
retrieving revision 1.1
diff -u -p -r1.1 rand.c
--- libgfortran/intrinsics/rand.c	12 Jun 2004 17:34:47 -0000	1.1
+++ libgfortran/intrinsics/rand.c	13 Jun 2004 20:01:27 -0000
@@ -77,17 +77,10 @@ prefix(irand) (GFC_INTEGER_4 *i)
 }
 
 
-/*  Return a REAL in the range [0,1).  Cast to double to use the full
-    range of pseudo-random numbers returned by irand().  */
+/*  Return a random REAL in the range [0,1).  */
 
 GFC_REAL_4
 prefix(rand) (GFC_INTEGER_4 *i)
 {
-  GFC_REAL_4 val;
-
-  do 
-    val = (GFC_REAL_4)((double)(prefix(irand) (i) - 1) / (double) 
GFC_RAND_M1);
-  while (val == 1.0);
-
-  return val;
+  return normalize_r4_i4 (i - 1, GFC_RAND_M1);
 }
Index: libgfortran/intrinsics/random.c
===================================================================
RCS file: /var/cvsroot/gcc-cvs/gcc/libgfortran/intrinsics/random.c,v
retrieving revision 1.8
diff -u -p -r1.8 random.c
--- libgfortran/intrinsics/random.c	13 Jun 2004 18:25:53 -0000	1.8
+++ libgfortran/intrinsics/random.c	13 Jun 2004 20:24:51 -0000
@@ -458,16 +458,11 @@ 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);
-      /* Burn a random number, so the REAL*4 and REAL*8 functions
-         produce similar sequences of random numbers.  */
-      kiss = kiss_random_kernel ();
-    }
-  while (*x == 1.0);
-
+  kiss = kiss_random_kernel ();
+  /* Burn a random number, so the REAL*4 and REAL*8 functions
+     produce similar sequences of random numbers.  */
+  kiss_random_kernel ();
+  *x = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
 }
 
 /*  This function produces a REAL(8) value from the uniform distribution
@@ -479,14 +474,9 @@ 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 == 1.0);
-
+  kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
+  kiss += kiss_random_kernel ();
+  *x = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
 }
 
 /*  This function fills a REAL(4) array with values from the uniform
Index: libgfortran/runtime/normalize.c
===================================================================
RCS file: libgfortran/runtime/normalize.c
diff -N libgfortran/runtime/normalize.c
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ libgfortran/runtime/normalize.c	13 Jun 2004 22:42:58 -0000
@@ -0,0 +1,111 @@
+/* Nelper routines to convert from integer to real.
+   Copyright 2004 Free Software Foundation, Inc.
+   Contributed by Paul Brook.
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+#include <math.h>
+#include "libgfortran.h"
+
+/* These routines can be sensitive to excess precision, so should really be
+   compiled with -ffloat-store.  */
+
+/* Return the largest value less than one representable in a REAL*4.  */
+
+static inline GFC_REAL_4
+almostone_r4 ()
+{
+#ifdef HAVE_NEXTAFTERF
+  return nextafterf (1.0f, 0.0f);
+#else
+  /* The volatile is a hack to prevent excess precision on x86.  */
+  static volatile GFC_REAL_4 val = 0.0f;
+  GFC_REAL_4 x;
+
+  if (val != 0.0f)
+    return val;
+
+  val = 0.9999f;
+  do
+    {
+      x = val;
+      val = (val + 1.0f) / 2.0f;
+    }
+  while (val > x && val < 1.0f);
+  if (val == 1.0f)
+    val = x;
+  return val;
+#endif
+}
+
+
+/* Return the largest value less than one representable in a REAL*8.  */
+
+static inline GFC_REAL_8
+almostone_r8 ()
+{
+#ifdef HAVE_NEXTAFTER
+  return nextafter (1.0, 0.0);
+#else
+  static volatile GFC_REAL_8 val = 0.0;
+  GFC_REAL_8 x;
+
+  if (val != 0.0)
+    return val;
+
+  val = 0.9999;
+  do
+    {
+      x = val;
+      val = (val + 1.0) / 2.0;
+    }
+  while (val > x && val < 1.0);
+  if (val == 1.0)
+    val = x;
+  return val;
+#endif
+}
+
+
+/* Convert an unsigned integer in the range [0..x) into a
+   real the range [0..1).  */
+
+GFC_REAL_4
+normalize_r4_i4 (GFC_UINTEGER_4 i, GFC_UINTEGER_4 x)
+{
+  GFC_REAL_4 r;
+
+  r = (GFC_REAL_4) i / (GFC_REAL_4) x;
+  if (r == 1.0f)
+    r = almostone_r4 ();
+  return r;
+}
+
+
+/* Convert an unsigned integer in the range [0..x) into a
+   real the range [0..1).  */
+
+GFC_REAL_8
+normalize_r8_i8 (GFC_UINTEGER_8 i, GFC_UINTEGER_8 x)
+{
+  GFC_REAL_8 r;
+
+  r = (GFC_REAL_8) i / (GFC_REAL_8) x;
+  if (r == 1.0)
+    r = almostone_r8 ();
+  return r;
+}