This is the mail archive of the
fortran@gcc.gnu.org
mailing list for the GNU Fortran project.
Re: [gfortran] Implement irand(), rand(), and srand() for G77 compatibility
- From: Steve Kargl <sgk at troutmask dot apl dot washington dot edu>
- To: Richard Henderson <rth at twiddle dot net>
- Cc: Paul Brook <paul at codesourcery dot com>, fortran at gcc dot gnu dot org, gcc-patches at gcc dot gnu dot org
- Date: Mon, 31 May 2004 18:13:09 -0700
- Subject: Re: [gfortran] Implement irand(), rand(), and srand() for G77 compatibility
- References: <20040530185518.GA61701@troutmask.apl.washington.edu> <200405302310.20676.paul@codesourcery.com> <20040530234527.GA63095@troutmask.apl.washington.edu> <200405310050.45040.paul@codesourcery.com> <20040531012128.GA63448@troutmask.apl.washington.edu> <20040531215531.GA2298@twiddle.net>
On Mon, May 31, 2004 at 02:55:31PM -0700, Richard Henderson wrote:
> On Sun, May 30, 2004 at 06:21:28PM -0700, Steve Kargl wrote:
>>> You're forgetting that real*4 has significantly less than 32 bits
>>> precision, so is unable to accurately represent integers of this
>>> magnitude. After rounding the result will in the range [0-1]
>>
>> Thanks for the wack with the clue stick.
>>
>> The simple fix is
>>
>> return (GFC_REAL_4) ((double)(prefix(irand)(i)-1) / (double)GFC_RAND_M1);
>
> Except if you *really* don't want 1.0 as a result, this won't work either,
> because the double->float conversion will round. I'm pretty sure what you
> want is truncated rounding, but I'm uncertain how to give you that
> portably.
Yes, I want the range [0,1), ie., inclusive of 0 and exclusive of 1.
> Probably the easiest portable fix is
>
> GFC_REAL_4 r = (prefix(irand)(i) - 1) / (GFC_REAL_4)GFC_RAND_M1;
> if (r == 1)
> r = nextafterf (1, 0);
> return r;
>
> mod whatever markup you use to select nextafter{,f,l} based on the
> underlying type of GFC_REAL_4.
The random_number intrinsic uses something similar to
do {
r = (GFC_REAL_4) ( (double) (prefix(irand)(i) - 1)
/ (double) GFC_RAND_M1);
} while (r == 1.0)
which should handle the double->float rounding problem.
--
Steve