Web pdp-10.trailing-edge.com

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
STD	XWAC2,-3(XPDP)
STD	XWAC4,-1(XPDP)
ARTEST	2,1
RDERR	1,DISCRETE: wrong number of subscripts
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
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
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]
EXEC	ALOG.
MOVN	XWAC1,X0
IF
CAMN	XWAC2,[-1.0]
GOTO	FALSE
THEN
RAND	3
MAKEREAL
ST	X0,.YFARG	;[205]
EXEC	ALOG.
FMPR	XWAC2,X0
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
L	XWAC3,AR(LOW,1)
SUB	XWAC3,AR(UPP,1)
SUBI	XWAC3,1
HRL	XWAC2,XWAC3
LI	X1,0
LOOP
AS
AOBJN	XWAC2,TRUE
SA
FMPR	X0,X1
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
HRRE	XWAC7,AR(LOW,2)
HRRE	XWAC6,AR(UPP,2)
HRLI	XWAC5,XWAC6
LOOP
L	XWAC10,XWAC6
AS
ASH	XWAC6,-1
CAMN	XWAC6,XWAC7
GOTO	FALSE
CAMG	XWAC3,@XWAC5
GOTO	TRUE
L	XWAC7,XWAC6
GOTO	TRUE+2
SA
CAMLE	XWAC3,@XWAC5
IF	;[232] Still inside array b
CAMLE	XWAC6,AR(UPP,2)
GOTO	FALSE
THEN	;Could be next interval
CAMLE	XWAC3,@XWAC5
FI
SUB	XWAC6,AR(LOW,2)
TLZ	XWAC6,-1
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
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
HRLI	X1,XWAC3
LOOP
L	XWAC3,XWAC4
AS
ASH	XWAC5,-1	;j:=j//2
JUMPE	XWAC5,FALSE
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
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))
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]
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
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
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
AS
CAML	X16,[1.0]
GOTO	TRUE
SA
ST	X16,(XPDP)
ST	X16,.YFARG		;[205]
EXEC	ALOG.
MOVN	X0,X0
FSC	X0,1
FDVR	X0,(XPDP)
ST	X0,.YFARG		;[205]
EXEC	SQRT.
SUB	XPDP,[2,,2]
FMPR	X0,1(XPDP)
FMPR	X0,ARG2
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
CAMG	X0,ARG2
RDERR	14,RANDINT or UNIFORM: too long interval
FI
;end of [110]
L	X0,ARG2
SUB	X0,ARG1
MUL	X0,X1
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
CAMGE	X1,ARG2
RDERR	14,RANDINT or UNIFORM: too long interval
FI
;end of [110]
L	X1,ARG2
FSBR	X1,ARG1
FMPR	X1,X0
```