[gfortran] Implement irand(), rand(), and srand() for G77 compatibility

[Steve Kargl]

On Mon, May 31, 2004 at 12:50:41AM +0100, Paul Brook wrote:
> > irand(i) returns a value in the range [1,GFC_RAND_M-1].
> > irand(i) - 1 is in the range [0,GFC_RAND_M-2].  So, the
> > division above is in the range [0,1).
> 
> 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);

I'm not sure if this works on amd64, alpha, sparc64, etc.

A more complicate fix is 

#define PRECS	32

GFC_REAL_4 div32 (GFC_INTEGER_4 num, GFC_INTEGER_4 dem)
{
  GFC_REAL_4 val;
  mpf_t n, d;
  mpf_init2 (mpf_t n, PRECS);
  mpf_init2 (mpf_t d, PRECS);
  mpf_set_ui (n, num);
  mpf_set_ui (d, dem);
  mpf_div (n, n, d);
  val = (GFC_REAL_4) mpf_get_d(n);
  mpf_clear (n);
  mpf_clear (d);
  return val;
}

/* Return a REAL in the range [0,1).  */
GFC_REAL_4 prefix(rand) (GFC_INTEGER_4 *i)
{
  GFC_REAL_4 val;
  val = div32(prefix(irand)(i) - 1, GFC_RAND_M1); 
  return val;
}

For rand(), I think the simple fix is in order.  Also, it seems
that libgfortran does not use gmp.

-- 
Steve