82 lines
3.8 KiB
Text
82 lines
3.8 KiB
Text
* Monte Carlo methods 08/03/2017
|
|
MONTECAR CSECT
|
|
USING MONTECAR,R13 base register
|
|
B 72(R15) skip savearea
|
|
DC 17F'0' savearea
|
|
STM R14,R12,12(R13) save previous context
|
|
ST R13,4(R15) link backward
|
|
ST R15,8(R13) link forward
|
|
LR R13,R15 set addressability
|
|
LA R8,1000 isamples=1000
|
|
LA R6,4 i=4
|
|
DO WHILE=(C,R6,LE,=F'7') do i=4 to 7
|
|
MH R8,=H'10' isamples=isamples*10
|
|
ZAP HITS,=P'0' hits=0
|
|
LA R7,1 j=1
|
|
DO WHILE=(CR,R7,LE,R8) do j=1 to isamples
|
|
BAL R14,RNDPK call random
|
|
ZAP X,RND x=rnd
|
|
BAL R14,RNDPK call random
|
|
ZAP Y,RND y=rnd
|
|
ZAP WP,X x
|
|
MP WP,X x**2
|
|
DP WP,ONE ~
|
|
ZAP XX,WP(8) x**2 normalized
|
|
ZAP WP,Y y
|
|
MP WP,Y y**2
|
|
DP WP,ONE ~
|
|
ZAP YY,WP(8) y**2 normalized
|
|
AP XX,YY xx=x**2+y**2
|
|
IF CP,XX,LT,ONE THEN if x**2+y**2<1 then
|
|
AP HITS,=P'1' hits=hits+1
|
|
ENDIF , endif
|
|
LA R7,1(R7) j++
|
|
ENDDO , enddo j
|
|
CVD R8,PSAMPLES psamples=isamples
|
|
ZAP WP,=P'4' 4
|
|
MP WP,ONE ~
|
|
MP WP,HITS *hits
|
|
DP WP,PSAMPLES /psamples
|
|
ZAP MCPI,WP(8) mcpi=4*hits/psamples
|
|
XDECO R6,WC edit i
|
|
MVC PG+4(1),WC+11 output i
|
|
MVC WC,MASK load mask
|
|
ED WC,PSAMPLES edit psamples
|
|
MVC PG+6(8),WC+8 output psamples
|
|
UNPK WC,MCPI unpack mcpi
|
|
OI WC+15,X'F0' zap sign
|
|
MVC PG+31(1),WC+6 output mcpi
|
|
MVC PG+33(6),WC+7 output mcpi decimals
|
|
XPRNT PG,L'PG print buffer
|
|
LA R6,1(R6) i++
|
|
ENDDO , enddo i
|
|
L R13,4(0,R13) restore previous savearea pointer
|
|
LM R14,R12,12(R13) restore previous context
|
|
XR R15,R15 rc=0
|
|
BR R14 exit
|
|
RNDPK EQU * ---- random number generator
|
|
ZAP WP,RNDSEED w=seed
|
|
MP WP,RNDCNSTA w*=cnsta
|
|
AP WP,RNDCNSTB w+=cnstb
|
|
MVC RNDSEED,WP+8 seed=w mod 10**15
|
|
MVC RND,=PL8'0' 0<=rnd<1
|
|
MVC RND+3(5),RNDSEED+3 return rnd
|
|
BR R14 ---- return
|
|
PSAMPLES DS 0D,PL8 F(15,0)
|
|
RNDSEED DC PL8'613058151221121' linear congruential constant
|
|
RNDCNSTA DC PL8'944021285986747' "
|
|
RNDCNSTB DC PL8'852529586767995' "
|
|
RND DS PL8 fixed(15,9)
|
|
ONE DC PL8'1.000000000' 1 fixed(15,9)
|
|
HITS DS PL8 fixed(15,0)
|
|
X DS PL8 fixed(15,9)
|
|
Y DS PL8 fixed(15,9)
|
|
MCPI DS PL8 fixed(15,9)
|
|
XX DS PL8 fixed(15,9)
|
|
YY DS PL8 fixed(15,9)
|
|
PG DC CL80'10**x xxxxxxxx samples give Pi=x.xxxxxx' buffer
|
|
MASK DC X'40202020202020202020202020202120' mask CL16 15num
|
|
WC DS PL16 character 16
|
|
WP DS PL16 packed decimal 16
|
|
YREGS
|
|
END MONTECAR
|