Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-03 - decus/20-0078/rts/rd.mac
There is 1 other file named rd.mac in the archive. Click here to see a list.
	SUBTTL	DECLARATIONS

COMMENT \
AUTHOR:		CLAES WIHLBORG
VERSION:	4	[143,205,232,243]
FUNCTION:	RANDOM DRAWING ROUTINES
CONTENTS:	.RDDI	DISCRETE
		.RDDR	DRAW
		.RDER	ERLANG
		.RDHI	HISTD
		.RDHO	HISTO (UTILITY)
		.RDLI	LINEAR
		.RDNE	NEGEXP
		.RDNO	NORMAL
		.RDPO	POISSON
		.RDRA	RANDINT
		.RDUN	UNIFORM

The algorithms use the fundamental function XI(u), which implements a
basic drawing giving a pseudo-random number in the interval (0,1)
(floating point). UI(u) is the discrete basic drawing which yields a
number in the interval (0,2^35) (integer, implemented by the RAND macro).
XI(u)=UI(u)/2^35.
u is the random seed, which is updated by UI or XI.
Essentially, next u:=u * 5^15 MOD 2^35.
\

	SEARCH	SIMMAC
	SALL
	SUBTTL	DISCRETE (.RDDI)

	SEARCH	SIMRPA
	RDINIT	DI

;INTEGER PROCEDURE DISCRETE(A,U);NAME U;ARRAY A;INTEGER U;
;	BEGIN INTEGER I,J;REAL X;
;		X:=XI(U); I:=LSB; J:=USB-I;
;		FOR J:=J//2 WHILE J NEQ 0 DO
;			IF A(I+J) LEQ X THEN I:=I+J;
;		FOR I:=I,I+1 WHILE I LEQ USB DO
;			IF A(I) GRT X THEN GOTO L;
;	L:	DISCRETE:=I
;	END;

.RDDI:	PROC
	RAND	2
	MAKEREAL
	EXCH	XWAC1,ARG1
	ADD	XPDP,[4,,4]
	STD	XWAC2,-3(XPDP)
	STD	XWAC4,-1(XPDP)
	ARTEST	2,1
	RDERR	1,DISCRETE: wrong number of subscripts
	LF	XWAC2,ZARBAD(XWAC1)
	L	XWAC4,AR(LOW,1)
	L	XWAC5,AR(UPP,1)
	L	XWAC3,AR(UPP,1)
	SUB	XWAC5,XWAC4
	HRLI	XWAC2,XWAC1
	LOOP
		L	XWAC1,XWAC4
	AS
		ASH	XWAC5,-1
		JUMPE	XWAC5,FALSE
		ADD	XWAC1,XWAC5
		CAMGE	X0,@XWAC2
		GOTO	TRUE
		L	XWAC4,XWAC1
		GOTO	TRUE+1
	SA
	LOOP	AS
		CAMGE	X0,@XWAC2
		GOTO	FALSE
		CAMGE	XWAC1,XWAC3
		AOJA	XWAC1,TRUE
		ADDI	XWAC1,1
	SA
	LD	XWAC4,-1(XPDP)
	LD	XWAC2,-3(XPDP)
	SUB	XPDP,[4,,4]
	EXCH	XWAC1,RESULT
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	DRAW (.RDDR)

	SEARCH	SIMRPA
	RDINIT	DR

;BOOLEAN PROCEDURE DRAW(A,U);NAME U;REAL A;INTEGER U;
;	DRAW:=XI(U) LESS A;

.RDDR:	PROC
	RAND	2
	MAKEREAL
;[143]	CAML	X0,ARG1
	SKIPLE	ARG1	;[143]
	CAMLE	X0,ARG1	;[143]
	TDZA	X1,X1
	SETO	X1,
	ST	X1,RESULT
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	ERLANG (.RDER)

	SEARCH	SIMRPA
	RDINIT	ER
	EXTERN	ALOG.

;REAL PROCEDURE ERLANG(A,B,U);NAME U;REAL A,B;INTEGER U;
;	BEGIN	REAL P,AB;
;		AB:=A*B;
;		P:=1;
;		FOR B:=B-1 WHILE B GEQ 0 DO
;			P:=P*XI(U)
;		P:=LN(P);
;		IF B NEQ -1 THEN
;			P:=P-B*LN(XI(U));
;		ERLANG:=-P/AB
;	END;

.RDER:	PROC
	EXCH	XWAC1,ARG1
	EXCH	XWAC2,ARG2
	CAIG	XWAC1,0
	RDERR	2,ERLANG: first parameter not positive
	CAIG	XWAC2,0
	RDERR	3,ERLANG: second parameter not positive
	FMPR	XWAC1,XWAC2
	PUSH	XPDP,XWAC1
	MOVSI	XWAC1,(1.0)
	WHILE
		FSBRI	XWAC2,(1.0)
		JUMPL	XWAC2,FALSE
	DO
		RAND	3
		MAKEREAL
		FMPR	XWAC1,X0
	OD
	ST	XWAC1,.YFARG		;[205]
	LI	16,.YFADR		;[205]
	EXEC	ALOG.
	MOVN	XWAC1,X0
	IF
		CAMN	XWAC2,[-1.0]
		GOTO	FALSE
	THEN
		RAND	3
		MAKEREAL
		ST	X0,.YFARG	;[205]
		LI	X16,.YFADR	;[205]
		EXEC	ALOG.
		FMPR	XWAC2,X0
		FADR	XWAC1,XWAC2
	FI
	FDVR	XWAC1,(XPDP)
	SUB	XPDP,[1,,1]
	L	XWAC2,ARG2
	EXCH	XWAC1,RESULT
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	HISTD (.RDHI)

	SEARCH	SIMRPA
	RDINIT	HI

;INTEGER PROCEDURE HISTD(A,U);NAME U;ARRAY A;INTEGER U;
;	BEGIN REAL S;INTEGER T;
;		FOR T:=LSB STEP 1 UNTIL USB DO
;			S:=S+A(T);
;		S:=S*XI(U);
;		T:=LSB;
;	L:	S:=S-A(T);
;		IF S GT 0 THEN BEGIN T:=T+1;GOTO L END;
;		HISTD:=T
;	END;

.RDHI:	PROC
	RAND	2
	MAKEREAL
	EXCH	XWAC1,ARG1
	STACK	XWAC2
	STACK	XWAC3
	ARTEST	2,1
	RDERR	4,HISTD: wrong number of subscripts
	LF	XWAC2,ZARBAD(XWAC1)
	L	XWAC3,AR(LOW,1)
	ADD	XWAC2,XWAC3
	SUB	XWAC3,AR(UPP,1)
	SUBI	XWAC3,1
	HRL	XWAC2,XWAC3
	LI	X1,0
	LOOP
		FADR	X1,(XWAC2)
	AS
		AOBJN	XWAC2,TRUE
	SA
	FMPR	X0,X1
	LF	X1,ZARBAD(XWAC1)
	HRLI	X1,XWAC1
	L	XWAC1,AR(LOW,1)
	LOOP
		FSBR	X0,@X1
	AS
		CAILE	X0,0
		AOJA	XWAC1,TRUE
	SA
	UNSTK	XWAC3
	UNSTK	XWAC2
	EXCH	XWAC1,RESULT
	RETURN
	EPROC
	LIT
	PRGEND
	SEARCH	SIMRPA
	RDINIT	HO

Comment \

PROCEDURE Histo(a,b,c,d); ARRAY a,b; REAL c,d;
BEGIN
	INTEGER i,j,k;
	i:=USB(b); j:=LSB(b);
	FOR k:=(i+j)//2 WHILE k NE i DO
	    IF b(k) < c THEN i:=k ELSE j:=k;
	IF b(k) < c THEN k:=k+1;
	IF k <= USB(b) THEN
	BEGIN IF b(k) < c THEN k:=k+1; END;
	a(k+LSB(a)-LSB(b)):=a(k+LSB(a)-LSB(b)) + d;
END
\

.RDHO:	PROC
	ARTEST	5,1
	RDERR	5,HISTO: wrong number of subscripts
	ARTEST	5,2
	RDERR	5,HISTO: wrong number of subscripts
	L	XWAC5,AR(UPP,1)
	SUB	XWAC5,AR(LOW,1)
	SUBI	XWAC5,1
	L	XWAC6,AR(UPP,2)
	SUB	XWAC6,AR(LOW,2)
	CAME	XWAC5,XWAC6
	RDERR	6,HISTO: incompatible array parameters
	LF	XWAC5,ZARBAD(XWAC2)
	HRRE	XWAC7,AR(LOW,2)
	HRRE	XWAC6,AR(UPP,2)
	HRLI	XWAC5,XWAC6
	LOOP
		L	XWAC10,XWAC6
		ADD	XWAC6,XWAC7
	AS
		ASH	XWAC6,-1
		CAMN	XWAC6,XWAC7
		GOTO	FALSE
		CAMG	XWAC3,@XWAC5
		GOTO	TRUE
		L	XWAC7,XWAC6
		ADD	XWAC6,XWAC10
		GOTO	TRUE+2
	SA
	CAMLE	XWAC3,@XWAC5
	ADDI	XWAC6,1
	IF	;[232] Still inside array b
		CAMLE	XWAC6,AR(UPP,2)
		GOTO	FALSE
	THEN	;Could be next interval
		CAMLE	XWAC3,@XWAC5
		ADDI	XWAC6,1
	FI
	SUB	XWAC6,AR(LOW,2)
	TLZ	XWAC6,-1
	ADD	XWAC6,AR(LOW,1)
	ADD	XWAC6,OFFSET(ZARBAD)(XWAC1)
	FADRM	XWAC4,(XWAC6)
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	LINEAR (.RDLI)

	SEARCH	SIMRPA
	RDINIT	LI

;REAL PROCEDURE LINEAR(A,B,U);NAME U;ARRAY A,B;INTEGER U;
;	BEGIN INTEGER I,J;REAL X;
;		X:=XI(U);I:=LBA;J:=UBA-I;
;		FOR J:=J//2 WHILE J NEQ 0 DO
;			IF A(I+J) LESS X THEN I:=I+J;
;		FOR I:=I+1 WHILE A(I) LESS X AND I LE UBA DO;
;		IF A(I)-A(I-1) LE 0 THEN LINEAR := B(I-1) ELSE
;		LINEAR:=B(I-1)+(X-A(I-1))*(B(I)-B(I-1))/(A(I)-A(I-1))
;	END;

.RDLI:	PROC
	RAND	3
	MAKEREAL
	EXCH	XWAC1,ARG1
	EXCH	XWAC2,ARG2
	ADD	XPDP,[4,,4]
	STD	XWAC3,-3(XPDP)
	STD	XWAC5,-1(XPDP)
	ARTEST	3,1
	RDERR	7,LINEAR: wrong number of subscripts
	ARTEST	3,2
	RDERR	7,LINEAR: wrong number of subscripts
	L	XWAC4,AR(LOW,1)
	CAME	XWAC4,AR(LOW,2)
	RDERR	10,LINEAR: array bounds do not match
	L	XWAC5,AR(UPP,1)
	CAME	XWAC5,AR(UPP,2)
	RDERR	8,LINEAR: array bounds do not match
	SUB	XWAC5,XWAC4
	LF	X1,ZARBAD(XWAC1)
	HRLI	X1,XWAC3
	LOOP
		L	XWAC3,XWAC4
	AS
		ASH	XWAC5,-1	;j:=j//2
		JUMPE	XWAC5,FALSE
		ADD	XWAC3,XWAC5	;i+j
		CAMG	X0,@X1
		GOTO	TRUE
		L	XWAC4,XWAC3
		GOTO	TRUE+1
	SA
	LOOP	AS
		CAMG	X0,@X1
		GOTO	FALSE
		CAML	XWAC3,AR(UPP,1)
		GOTO	FALSE
		AOJA	XWAC3,TRUE
	SA
	L	XWAC4,XWAC3
	ADD	XWAC4,OFFSET(ZARBAD)(XWAC2)
	L	XWAC5,@X1	;A(i)
	L	XWAC6,(XWAC4)	;B(i)
	SUBI	XWAC3,1		;i-1 for A
	SUBI	XWAC4,1		;i-1 for B
	FSBR	X0,@X1		;X-A(i-1)
	FSBR	XWAC5,@X1	;A(i)-A(i-1)
	IF	;A(I)-A(I-1) > 0
		JUMPLE	XWAC5,FALSE
	THEN	FSBR	XWAC6,(XWAC4)	;B(i)-B(i-1)
		FMPR	X0,XWAC6	; * (X-A(i-1))
		FDVR	X0,XWAC5	; /(A(i)-A(i-1))
		FADR	X0,(XWAC4)	; + B(i-1)
	ELSE	;Result is B(i-1)
		L	X0,(XWAC4)
	FI
	LD	XWAC5,-1(XPDP)
	LD	XWAC3,-3(XPDP)
	SUB	XPDP,[4,,4]
	L	XWAC2,ARG2
	L	XWAC1,ARG1
	ST	X0,RESULT
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	NEGEXP (.RDNE)

	SEARCH	SIMRPA
	RDINIT	NE
	EXTERN	ALOG.

;REAL PROCEDURE NEGEXP(A,U);NAME U;REAL A;INTEGER U;
;	NEGEXP:=-LN(XI(U))/A;

.RDNE:	PROC
	SKIPG	ARG1
	RDERR	11,NEGEXP: first argument not positive
	RAND	2
	MAKEREAL
	ST	X0,.YFARG		;[205]
	LI	X16,.YFADR		;[205]
	EXEC	ALOG.
	MOVN	X0,X0
	FDVRM	X0,ARG1
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	POISSON (.RDPO)

	SEARCH	SIMRPA
	RDINIT	PO
	EXTERN	SQRT.,EXP.,.RDNO

;INTEGER PROCEDURE POISSON(A,U);NAME U;REAL A;INTEGER U;
;	BEGIN INTEGER T;REAL R,R1;
;
;		IF A<=25 THEN BEGIN
;			R1:=EXP(-A);R:=1;
;			L: R:=R*XI(U); IF R GEQ R1 THEN
;			BEGIN T:=T+1;GOTO L END;
;			END
;		ELSE
;			T:=MAX(0,ENTIER(NORMAL(A,SQRT(A),U)+0.5));
;		POISSON:=T
;	END

.RDPO:	PROC
	MOVN	X0,ARG1
	LI	X16,.YFADR		;[205]
	IF
		CAMGE	X0,[-25.0]
		GOTO	FALSE
	THEN
		ST	X0,.YFARG	;[205]
		EXEC	EXP.
		STACK	X0
		ST	XWAC1,ARG1
		LI	XWAC1,0
		MOVSI	X16,(1.0)
		LOOP
			RAND	2
			MAKEREAL
			FMPR	X16,X0
		AS
			CAML	X16,(XPDP)
			AOJA	XWAC1,TRUE
		SA
		UNSTK	X0
	ELSE
		MOVNM	X0,.YFARG	;[205]
		EXEC	SQRT.
		EXCH	XWAC1,ARG1
		EXCH	XWAC3,ARG2
		STACK	XWAC2
		STACK	XTAC
		LI	XTAC,XWAC1
		ST	X0,ARG2
		EXEC	.RDNO
		FIXR	XWAC1,XWAC1
		SKIPGE	XWAC1
		LI	XWAC1,0
		UNSTK	XTAC
		UNSTK	XWAC2
		EXCH	XWAC3,ARG2
	FI
	EXCH	XWAC1,RESULT
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	NORMAL (.RDNO)

	SEARCH	SIMRPA
	RDINIT	NO
	EXTERN	ALOG.,SQRT.

Comment \ Algorithm:

REAL PROCEDURE Normal(a,b,u);NAME u;REAL a,b;INTEGER u;
	BEGIN REAL x,s;
	L:	x:=2*XI(u)-1;
		s:=(2*XI(u)-1)**2+x*x;
		IF s >= 1 THEN GOTO L;
		Normal:=a+b*x*Sqrt(-2*Ln(s)/s)
	END;
\

.RDNO: PROC
	SKIPG	ARG2
	RDERR	12,NORMAL:  standard deviation not positive
	ADD	XPDP,[2,,2]
	LOOP
		RAND	3
		MAKEREAL(2)
		FSBRI	X0,(1.0)
		ST	X0,-1(XPDP)
		FMPR	X0,X0
		L	X16,X0
		RAND	3
		MAKEREAL (2)		;[243]
		FSBRI	X0,(1.0)	;[243]
		FMPR	X0,X0
		FADR	X16,X0
	AS
		CAML	X16,[1.0]
		GOTO	TRUE
	SA
	ST	X16,(XPDP)
	ST	X16,.YFARG		;[205]
	LI	X16,.YFADR		;[205]
	EXEC	ALOG.
	MOVN	X0,X0
	FSC	X0,1
	FDVR	X0,(XPDP)
	ST	X0,.YFARG		;[205]
	LI	X16,.YFADR		;[205]
	EXEC	SQRT.
	SUB	XPDP,[2,,2]
	FMPR	X0,1(XPDP)
	FMPR	X0,ARG2
	FADRM	X0,ARG1
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	RANDINT (.RDRA)

	SEARCH	SIMRPA
	RDINIT	RA

;INTEGER PROCEDURE RANDINT(A,B,U);NAME U;INTEGER A,B,U;
;	RANDINT:=UI(U)*(B-A+1)//2**35+A;

.RDRA:	PROC
	RAND	3
	;[110] check interval
	L	X0,ARG1
	CAMLE	X0,ARG2
	RDERR	13,RANDINT or UNIFORM: interval negative
	IF
		;check abs value B-A<2**35-1
		JUMPG	X0,FALSE			;interval ok
	THEN
		;A is not positive, B might be
		;check length of interval
		ADD	X0,[377777,,-1]			;Add greatest number
		CAMG	X0,ARG2
		RDERR	14,RANDINT or UNIFORM: too long interval
	FI
	;end of [110]
	L	X0,ARG2
	SUB	X0,ARG1
	ADDI	X0,1
	MUL	X0,X1
	ADDM	X0,ARG1
	RETURN
	EPROC
	LIT
	PRGEND
	SUBTTL	UNIFORM (.RDUN)

	SEARCH	SIMRPA
	RDINIT	UN

;REAL PROCEDURE UNIFORM(A,B,U);NAME U;REAL A,B;INTEGER U;
;	UNIFORM:=XI(U)*(B-A)+A;

.RDUN:	PROC
	RAND	3
	MAKEREAL
	;[110] check interval
	L	X1,ARG1
	CAMLE	X1,ARG2
	RDERR	13,RANDINT or UNIFORM: interval negative
	IF
		;check abs B-A<=max real
		JUMPGE X1,FALSE			;interval ok
	THEN
		;A is negative, B might be
		;check length of interval
		FAD	X1,[377777,,-1]			;add greates number
		CAMGE	X1,ARG2
		RDERR	14,RANDINT or UNIFORM: too long interval
	FI
	;end of [110]
	L	X1,ARG2
	FSBR	X1,ARG1
	FMPR	X1,X0
	FADRM	X1,ARG1
	RETURN
	EPROC
	LIT
	END