.MCALL	.MODULE
.MODULE RAN$,VERSION=11,COMMENT=<RANDOM NUMBER SERVICE ROUTINE>,LIB=YES,IDENT=NO

; Copyright (c) 1998 by Mentec, Inc., Nashua, NH.
; All rights reserved
;
; This software is furnished under a license for use only on a
; single computer system and  may  be  copied  only  with  the
; inclusion of the above copyright notice.  This software,  or
; any other copies thereof, may not be provided  or  otherwise
; made available to any other person except for  use  on  such
; system and to one who agrees to these license  terms.  Title
; to  and  ownership of the software shall at all times remain
; in Mentec, Inc.
;
; The information in this document is subject to change without
; notice and should not be construed as a commitment by Digital
; Equipment Corporation, or Mentec, Inc.
;
; Digital and Mentec  assume  no  responsibility for the use or
; reliability of its software on equipment which is not supplied
; by Digital  or Mentec,  and listed  in the  Software  Product
; Description.
;
; AUTHOR:
;   R. SCHAEFER	24-OCT-78
;	DEVELOPED FROM THE VAX-RTL VERSION DUE TO:
;	P.F. CONKLIN, T.HASTINGS, D. MACLAREN
;
;+
; FUNCTION:
;	RANDOM NUMBER GENERATOR
;
; CALLING SEQUENCE:
;	MOV	#SEED,R0
;	JSR	PC,RAN$
;
;	RETURNS A PSEUDO-RANDOM NUMBER IN THE RANGE [0.0,1.0)
;
;	PRESERVES R5, DESTROYS ALL OTHER REGISTERS
;
; ALGORITHM:
;	THIS IS A GENERAL RANDOM NUMBER GENERATOR.  IT IS
;	OF THE MULTIPLICATIVE CONGRUENTIAL TYPE, AND HENCE
;	IS FAST, ALTHOUGH PRONE TO CERTAIN CLASSES OF
;	NON-RANDOM SEQUENCES.  THIS NON-RANDOM BEHAVIOR
;	TYPICALLY ARISES WHEN CONSIDERING TRIPLES OF NUMBERS
;	GENERATED BY THIS METHOD.
;
;	FOR MORE INFORMATION, CONSULT:
;	RANDOM NUMBER GENERATION (PP. 1192-1197)
;	BY G. MARSAGLIA IN
;	ENCYCLOPEDIA OF COMPUTER SCIENCE
;	EDITED BY ANTHONY RALSTON
;	PETROCELLI (NEW YORK 1976)
;
;METHOD:
;	SEED = 69069.*SEED + 1  (MOD 2**32)
;	69069. =  6715,1 AS 2 OCTAL WORDS
;
;	THE MULTIPLY ALGORITHM IS THAT USED BY THE PDP-10 HARDWARE
;	MULTIPLY.  THE COMMENTING REFLECTS INTERNAL REGISTERS.
;	IT IS SUGGESTED THAT THE READER REFER TO THE -10
;	ASSEMBLY LANGUAGE MANUAL IN READING THIS CODE.
;-
.IIF NDF EIS$I, EIS$I=0
.IF EQ EIS$I
	.MCALL	SOB		;RH001
.ENDC; EQ EIS$I

	.PSECT	SYS$I,I		;RH001

RAN$::
	MOV	#6715,R2	;INIT LOW ORDER MQ
	MOV	#1,R1		;INIT HIGH ORDER MQ
	MOV	R5,-(SP)	;SAVE R5
	MOV	#32.,R5		;COUNT
	CLR	R3		;CLEAR AR AND C-BIT
	CLR	R4
10$:	MOV	R2,-(SP)	;TEST IF BIT SHIFTED INTO MQ<0,1>
	ADC	@SP		;IS SAME AS BIT SHIFTED OUT (C-BIT)
	ROR	(SP)+
	BCC	30$
	BIT	#1,R2		;TEST MQ<0,1>
	BEQ	20$
	SUB	(R0),R4		;SUB BR FROM AR
	SBC	R3
	SUB	2(R0),R3
	BR	30$
20$:	ADD	(R0),R4		;ADD BR TO AR
	ADC	R3
	ADD	2(R0),R3
30$:	ASR	R3		;DOUBLE SHIFT AR AND MQ
	ROR	R4
	ROR	R1
	ROR	R2		;BIT SHIFTED OUT INTO C-BIT
	SOB	R5,10$		;LOOP
	MOV	(SP)+,R5	;RESTORE R5
	ADD	#1,R2		;INCR TO PROTECT AGAINST STRANGE SEEDS
	ADC	R1
	MOV	R2,(R0)+	;STORE NEW LOW ORDER SEED
	MOV	R1,(R0)		;STORE NEW HIGH ORDER SEED

; NOW CONVERT TO FLOATING POINT
;
	CLRB	R2		;SAVE HIGH ORDER 24 RESULT BITS
	MOV	#200,R0		;GET INITIAL EXPONENT
NORM:	ASL	R2		;FLOAT RESULT
	ROL	R1
	BCS	EXP		;JUMP IF LEADING BIT FOUND
	DEC	R0		;COMPENSATE EXPONENT
	BNE	NORM

; NOW HAVE NORMALIZED, OR GOT 0.0
;
EXP:	SWAB	R0		;NOW MAKE INTO FLOATING POINT FORMAT
	SWAB	R1
	SWAB	R2
	BISB	R1,R0
	CLRB	R1
	BISB	R2,R1
	ROR	R0		;INSERT + SIGN
	ROR	R1
	RETURN

	.END