From:	SMTP%"RELAY-INFO-VAX@CRVAX.SRI.COM" 14-MAR-1994 08:19:53.22
To:	EVERHART
CC:	
Subj:	RE: Random numbers.

From: harvey@indyvax.iupui.edu (James Harvey)
X-Newsgroups: comp.os.vms
Subject: RE: Random numbers.
Message-Id: <1994Mar10.215924.4201@ivax>
Date: 10 Mar 94 21:59:24 -0500
Organization: Indiana University-Purdue University at Indianapolis
Lines: 218
To: Info-VAX@CRVAX.SRI.COM
X-Gateway-Source-Info: USENET

In article <0097B2E4.B0C69140.3914@nscvax.princeton.edu>, Mighty Firebreather <dragon@nscvax.princeton.edu> writes:
> 	dh6375@albnyvms.bitnet (Master Dan) writes:
>
>>Can I make psuedo-random numbers any other way then using f$time?
>>
>
> 	Yes.
>
> 	Hint:  HELP RTL MTH$ MTH$RANDOM.
>
> 	Additional hint: test it to make sure it's adequately random for
> your purposes.

At least the last time I checked about five years ago, this one was pretty
horrible.  As I recall it is worse than the standard UNIX one (see below).

> 	Extra additional hint: see Knuth, Donald E. "The Art of Computer
> Programming: Volume 2 - Semi-Numerical Algorithms" on the subject.  Damned
> few pseudo random nubmer generators produce "good" sequences of random
> numbers.  As an example, the VAX C random number generator was criticized
> here recently for alternately producing odd and even numbers.  This was
> explained as being "bug for bug" compatibility with Unix!

Yep.  The "standard" UNIX random number generator is worthless for anything
but games.  And I wouldn't even bet on that!

>                                                            Another
> generator has the interesting property that if successive "random" values
> are plotted as points in 3D, they will fall in distinct planes.  And so it
> goes. <sigh>

The moral of the story is, if you need a real good one, use two fairly
good generators (that hopefully use different techniques) and shuffle
(the shuffling technique is described in Knuth) or hire a mathemetician!

Here is an implementation in MACRO of a "Minimal Standard Pseudo-Random
Number Generator" that was described in an article in the October 1988
issue of CACM.  As I recall it was supposed to be better than your average
maximal-length linear congruential generator (it is not a maximal-length
generator, but generates a sequence that is nearly maximal length, and it
is very fast as well).  No warranties expressed or implied of course.  HTH.

=====8<==clip'n'save==8<==clip'n'save==8<==clip'n'save==8<==clip'n'save==8<=====
	.TITLE	RANDPACK - Minimal Standard Pseudo-Random Number Generator
	.IDENT	'V00-000'

;++
; FACILITY:
;
;	VAX/VMS Minimal Standard Pseudo-Random Number Generator
;
; ABSTRACT:
;
;	This package contains an efficient assembly-language implementation
;	of the "minimal standard" random number generator similar to the
;	one described in the article "Random: Number Generators: Good Ones
;	are Hard to Find" by Stephen K. Park and Keith W. Miller, CACM
;	vol. 31 no. 10 (Oct. 1988), pp 1192-1201.
;
; AUTHOR:
;
;	06-Feb-1989	James A. Harvey
;
;	This program has been granted to the public domain by the author.
;
; Revision History:
;
;	None (yet).
;

.SBTTL	Symbolic Constants and Impure Data

	MODULUS=2147483647		; The modulus, m.
	FMODULUS=^F2147483647.0		; Ditto, as F-floating number.
	MINSEED=1			; The minimum value for SEED.
	MAXSEED=2147483646		; The maximum value for SEED.
	MULTIPLIER=16807		; The multiplier, a.

.PSECT	RWDATA,LONG,NOSHR,WRT		; Impure data.

SEED:	.LONG	1			; The default seed.

TIME:	.BLKQ	1			; For calls to RANDOMIZE(0)



.SBTTL	RANDOMIZE - Initialize the Pseudo-Random Number Generator
;++
;
; FUNCTIONAL DESCRIPTION:
;
;	RANDOMIZE [seed]
;
;	Initializes the pseudo-random number generator with the seed
;	passed, if seed is in the range 1 to 2**31-2.  If the seed
;	argument is omitted or not in the valid range, the generator's
;	seed is initialized based on the time of day as returned by the
;	$GETTIM system service.
;
; CALLING SEQUENCE:
;
;	Called with CALLS or CALLG to RANDOMIZE.  The argument is
;	optional.
;
; INPUT PARAMETERS:
;
;	seed
;	VMS Usage:	longword_unsigned
;	type:		longword (unsigned)
;	access:		read only
;	mechanism:	by reference
;
;	Seed is the initial seed for the random number generator.  If
;	the seed passed is not in the valid range of acceptable seeds
;	(1 to 2**31-2) the generator's seed is initialized with a value
;	based on the current system time as returned by the $GETTIM
;	system service.
;
; RETURN VALUE:
;
;	VMS Usage:	longword_unsigned
;	type:		longword (unsigned)
;	access:		write only
;	mechanism:	by value
;
;	RANDOMIZE returns the seed used to initialize the random number
;	generator.  This is always the same as the seed argument passed
;	to RANDOMIZE if the argument was in the range 1 to 2**31-2.  If
;	the argument was omitted, or not in the required range, the
;	returned value is the seed actually used to initialize the random
;	number generator.
;
; SIDE EFFECTS:
;
;	Updates SEED in impure data area.
;	All registers except R0 and R1 are preserved.
;

.PSECT	CODE,LONG,SHR,NOWRT

	NARGS=0				; Offset to number of arguments.
	SEEDARG=4			; Offset to optional seed argument.

.ENTRY	RANDOMIZE,<0>
	TSTB	NARGS(AP)		; Any argument?
	BEQLU	10$			; No - randomize by time clock.
	MOVL	SEEDARG(AP),R0		; Move argument to R0.
	CMPL	R0,#MINSEED		; Is argument less than 1?
	BLSS	10$			; If so, it's out of range.
	CMPL	R0,#MAXSEED		; Is argument greater than 2**31-2?
	BLEQ	30$			; If not, use it.
10$:	$GETTIM	TIMADR=TIME		; Get the system time in TIME:TIME+7
	MOVL	TIME,R0			; Get high order word in R0.
	XORL2	TIME+4,R0		; Exclusive or in least sig. bits.
	BICL2	#^X80000000,R0		; Clear the sign bit.
	CMPL	R0,#MINSEED		; Is argument less than 1?
	BLSS	20$			; If so, it's out of range, use 1.
	CMPL	R0,#MAXSEED		; Is argument greater than 2**31-2?
	BLEQ	30$			; If not, use it.
20$:	MOVL	#1,R0			; Use 1 if $GETTIM result out of range.
30$:	MOVL	R0,SEED			; Store the initial seed.
	RET				; And return.



.SBTTL	RANDOM - Return Next Pseudo-Random Uniform Variate on 0..1
;++
;
; FUNCTIONAL DESCRIPTION:
;
;	RANDOM
;
;	This routine returns an F_floating pseudo-random number uniformly
;	distributed in the range 0..1.
;
; CALLING SEQUENCE:
;
;	Called with CALLS or CALLG to RANDOM.  No arguments are required.
;
; INPUT PARAMETERS:
;
;	None.
;
; RETURNS:
;
;	VMS Usage:	floating_point
;	type:		F_floating
;	access:		write only
;	mechanism:	by value
;
;	RANDOM routine returns an F_floating pseudo-random number uniformly
;	distributed in the range 0..1.
;
; SIDE EFFECTS:
;
;	Updates SEED in impure data area.
;	All registers except R0 and R1 are preserved.
;

.PSECT	CODE,LONG,SHR,NOWRT

.ENTRY	RANDOM,<0>
	EMUL	#MULTIPLIER,SEED,#0,R0	; R0 = a * z
	EDIV	#MODULUS,R0,SEED,SEED	; z = (a * z) mod m.
	CVTLF	SEED,R0			; R0 = float(z)
	DIVF2	#FMODULUS,R0		; R0 = float(z) / float(m).
	RET

.ENTRY	GETSEED,<0>
	MOVL	SEED,R0
	RET
.END
=====8<==clip'n'save==8<==clip'n'save==8<==clip'n'save==8<==clip'n'save==8<=====
--
James Harvey   harvey@iupui.edu   IUPUI OIT Technical Support/VMS/Unix/Networks
Disclaimer:  These are my own opinions.  I do not speak for Indiana University.
Mathematicians do it in theory.


