#include <math.h>

/*... function u16807d() returns a U(0,1) random number ...................
c..........................................................................
c     bruce schmeiser     january 1992.                                   .
c     a linear congruential pseudorandom number generator                 .
c       using constant 16807 and modulus (2**31)-1.                       .
c              iseed = iseed*16807 (mod 2^31 -1)                          .
c     for implementations that don't require double precision, see        .
c       s.k. park and k.w. miller, "random numbers generators: good       .
c         ones are hard to find," cacm, 31, 10 (October 1988), 1192-1201. .
c     in correct implementations, starting with a seed value of 1 will    .
c       result in a seed value of 1043618065 on call number 10001.        .
c     note:  in C, real values are by default double precision.
c..........................................................................
c     input:  iseed.   integer.                                           .
c                        chosen from [1,2147483646] on the first call.    .
c                        thereafter, the value returned from the last call.
c     output: iseed.   integer.                                           .
c                        to be used in the next call.                     .
c     output: u16807d.  real.                                             .
c                        a pseudorandom number in (0,1).                  .
c........................................................................*/
float u16807d( iseed )
int *iseed;
{
      *iseed = fmod( *iseed * 16807., 2147483647. );
      return ( *iseed / 2147483647. );
}