Trailing-Edge
-
PDP-10 Archives
-
decuslib10-05
-
43,50337/16/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