+---------------------------------+
|              RAND               |
+---------------------------------+
RAND([<expN>])

-----------------------------------
Returns uniformly distributed random number between 0 and 1.
Return value - Numeric
-----------------------------------

The RAND() function returns a random number between 0 and 1.  Used
without an argument, FoxPro automatically applies a seed of 100001 to
the function.  You can restore the default at any time with
RAND(100001).

Including an optional numeric expression <expN> allows you to specify
the seed value to use with the function.  Using the same positive value
for <expN> always returns the same sequence of random numbers.

If <expN> is negative, a seed value is created from the system clock.
To achieve maximum randomness, use RAND() initially with a negative
argument and subsequent RAND()'s without an argument.

+---------------------------------+
|            Examples             |
+---------------------------------+

* ---------------------------------
*       Seed random number generator
*       (do this only once per program)
* ---------------------------------
= RAND(-1)

* ---------------------------------
*       Print uniformly distributed
*   random number between 0 and 1
* ---------------------------------
? RAND()

* ---------------------------------
*       Returns a uniformly distributed
*       real number between a and b
* ---------------------------------
FUNCTION urand
PARAMETER a, b
RETURN (b-a)*rand()+a

* ---------------------------------
*       Returns a uniformly distributed
*       integer between i and j
* ---------------------------------
FUNCTION irand
PARAMETER i, j
RETURN int((j-i+1)*rand()+i)

* ---------------------------------
*       Returns a value which simulates
*       throwing a pair of dice.
*       (Uses 'FUNCTION irand' above.)
*
*       Note:  This is NOT the same as
*       'irand(12)' which has a quite
*       different probability
*       distribution.
* ---------------------------------
FUNCTION dice
RETURN irand(1,6)+irand(1,6)

* ---------------------------------
*       Returns a random n-digit number
*       (Uses 'FUNCTION irand' above)
* ---------------------------------
FUNCTION digits
PARAMETER n
x = 10
y = irand(0, x-1) + x
RETURN substr(str(y,n+1),2)

* ---------------------------------
* Function to return exponentially
* distributed values of mean 'a'
* ---------------------------------
FUNCTION edist
PARAMETER a
RETURN -a*log(rand())

* ---------------------------------
*       Function to return a normally
*       distributed value with mean 'm'
*       and standard deviation 'std'
* ---------------------------------
*       This function requires that
*       the subroutine 'x_p' also be
*       present.
*
*       The algorithm was taken from
*       the "Handbook of Mathematical
*       Functions,"National Bureau of
*       Standards, June 1964.  Equation
*       26.2.23, page 933.
*
*       Its error is < (4.5E-4 * std)
*
FUNCTION ndist
PARAMETER m, std
x = x_p(RAND())
RETURN x * std + m

FUNCTION x_p
PARAMETER p
IF p < 0  or p > 1
        RETURN 0
ENDIF
IF p < .5
        adjust = .T.
ELSE
        p = 1 - p
        adjust = .F.
ENDIF
t = SQRT(LOG(1/(p*p)))
x = t -((t * 0.010328 + 0.802853) ;
                * t + 2.515517) ;
        / (((t * 0.001308 + 0.189269) ;
                * t + 1.432788) * t + 1)
RETURN IIF(adjust, -x, x)

* ---------------------------------
*       Returns area under normal curve
*       with mean 'm' and standard
*       deviation 'std' from -. to x
* ---------------------------------
*       The algorithm was taken from
*       the "Handbook of Mathematical
*       Functions,"National Bureau of
*       Standards, June 1964.
*   Equation 26.2.17, page 932.
*       Its error is < (7.5E-8 * std)
*
*       This doesn't use the RAND()
*       function, it's just included
*       here for your convenience.
*
*       It requires the function 'p'
*       to also be present
*
FUNCTION normal
PARAMETERS x, m, std
RETURN p((x-m)/std)

FUNCTION p
PARAMETER x
IF x < 0
        adjust = .T.
        x = -x
ELSE
        adjust = .F.
ENDIF
t = 1 / (1 + 0.2316419 * x)
poly = ((((1.330274429 ;
                * t - 1.821255978) ;
                * t + 1.781477937) ;
                * t - 0.356563782) ;
                * t + 0.319381530) * t
res =  1 - 0.398942280401432678 ;
                        * EXP(-x*x/2) * poly