Monte Carlo methods: Difference between revisions
Content added Content deleted
(add RPL) |
(Added solution for EDSAC) |
||
Line 1,042: | Line 1,042: | ||
3.1407 |
3.1407 |
||
3.1413 |
3.1413 |
||
=={{header|EDSAC order code}}== |
|||
Because real numbers on EDSAC were restricted to the interval [-1,1), this solution estimates pi/10 instead of pi. With 100,000 trials the program would have taken about 3.5 hours on the original EDSAC. |
|||
<syntaxhighlight lang="edsac"> |
|||
[Monte Carlo solution for Rosetta Code.] |
|||
[EDSAC program, Initial Orders 2.] |
|||
[Arrange the storage] |
|||
T45K P56F [H parameter: library s/r P1 to print real number] |
|||
T46K P78F [N parameter: library s/r P7 to print integer] |
|||
T47K P210F [M parameter: main routine] |
|||
T48K P114F [& (delta) parameter: library s/r C6 (division)] |
|||
T49K P150F [L parameter: library subroutine R4 to read data] |
|||
T51K P172F [G parameter: generator for pseudo-random numbers] |
|||
[Library subroutine M3, runs at load time and is then overwritten. |
|||
Prints header; here the header sets teleprinter to figures.] |
|||
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF |
|||
*!!!!TRIALS!!EST!PI#X10@&.. |
|||
PK [after header, blank tape and PK (WWG, 1951, page 91)] |
|||
[================ Main routine ====================] |
|||
E25K TM GK |
|||
[Variables] |
|||
[0] PF PF [target count: print result when count = target] |
|||
[2] PF PF [count of points] |
|||
[4] PF PF [count of hits (point inside circle)] |
|||
[6] PF PF [x-coordinate - 1/2] |
|||
[Constants] |
|||
T8#Z PF T10#Z PF [clear sandwich bits in 35-bit constants] |
|||
T8Z [resume normal loading] |
|||
[8] PD PF [35-bit constant 1] |
|||
[10] L1229F Y819F [35-bit constant 2/5 (near enough)] |
|||
[12] IF [1/2] |
|||
[13] RF [1/4] |
|||
[14] #F [figures shift] |
|||
[15] MF [dot (decimal point) in figures mode] |
|||
[16] @F [carriage return] |
|||
[17] &F [line feed] |
|||
[18] !F [space] |
|||
[Enter with acc = 0] |
|||
[19] A19@ GL [read seed for LCG into 0D] |
|||
AD T4D [pass seed to LCG in 4D] |
|||
[23] A23@ GG [initialize LCG] |
|||
T2#@ T4#@ [zero trials and hits] |
|||
[Outer loop: round target counts] |
|||
[27] TF [clear acc] |
|||
[28] A28@ GL [read next target count into 0D] |
|||
SD [acc := -target] |
|||
E85@ [exit if target = 0] |
|||
T#@ [store negated target] |
|||
[Inner loop : round points in the square] |
|||
[33] TF T4D [pass LCG range = 0 to return random real in [0,1)] |
|||
[35] A35@ G1G [call LCG, 0D := random x] |
|||
AD S12@ T6#@ [store x - 1/2 over next call] |
|||
T4D |
|||
[41] A41@ G1G [call LCG, 0D := random y] |
|||
AD S12@ TD [store y - 1/2] |
|||
H6#@ V6#@ [acc := (x - 1/2)^2] |
|||
HD VD [acc := acc := (x - 1/2)^2 + (y - 1/2)^2] |
|||
S13@ [test for point inside circle, i.e. acc < 1/4] |
|||
E56@ [skip if not] |
|||
TF A4#@ A8#@ T4#@ [inc number of hits] |
|||
[56] TF A2#@ A8#@ U2#@ [inc number of trials] |
|||
A#@ [add negated target] |
|||
G33@ [if not reached target, loop back] |
|||
A2#@ TD [pass number of trials to print s/r] |
|||
[64] A64@ GN [print number of trials] |
|||
A4#@ TD A2#@ T4D [pass hits and trials to division s/r] |
|||
[70] A70@ G& [0D := hits/trials, estimated value of pi/4] |
|||
HD V10#@ TD [times 2/5; pass estimated pi/10 to print s/r] |
|||
O18@ O18@ O8@ O15@ [print ' 0.'] |
|||
[79] A79@ GH P5F [print estimated pi/10 to 5 decimals] |
|||
O16@ O17@ [print CR, LF] |
|||
E27@ [loop back for new target] |
|||
[85] O14@ [exit: print dummy character to flush printer buffer] |
|||
ZF [halt program] |
|||
[==================== Generator for pseudo-random numbers ===========] |
|||
[Linear congruential generator, same algorithm as Delphi 7 LCG. |
|||
38 locations] |
|||
E25K TG |
|||
GK G10@ G15@ T2#Z PF T2Z I514D P257F T4#Z PF T4Z PD PF T6#Z PF T6Z PF RF A6#@ S4#@ T6#@ E25F E8Z PF T8Z PF PF A3F T14@ A4D T8#@ ZF A3F T37@ H2#@ V8#@ L512F L512F L1024F A4#@ T8#@ H6#@ C8#@ T8#@ S4D G32@ TD A8#@ E35@ H4D TD V8#@ L1F TD ZF |
|||
[==================== LIBRARY SUBROUTINES ============================] |
|||
[D6: Division, accurate, fast. |
|||
36 locations, workspace 6D and 8D. |
|||
0D := 0D/4D, where 4D <> 0, -1.] |
|||
E25K T& GK |
|||
GKA3FT34@S4DE13@T4DSDTDE2@T4DADLDTDA4DLDE8@RDU4DLDA35@ |
|||
T6DE25@U8DN8DA6DT6DH6DS6DN4DA4DYFG21@SDVDTDEFW1526D |
|||
[R4: Input of one signed integer at runtime. |
|||
22 storage locations; working positions 4, 5, and 6.] |
|||
E25K TL |
|||
GKA3FT21@T4DH6@E11@P5DJFT6FVDL4FA4DTDI4FA4FS5@G7@S5@G20@SDTDT6FEF |
|||
[P1: Prints non-negative fraction in 0D, without '0.'] |
|||
E25K TH |
|||
GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F |
|||
[P7, prints long strictly positive integer; |
|||
10 characters, right justified, padded left with spaces. |
|||
Even address; 35 storage locations; working position 4D.] |
|||
E25K TN |
|||
GKA3FT26@H28#@NDYFLDT4DS27@TFH8@S8@T1FV4DAFG31@SFLDUFOFFFSF |
|||
L4FT4DA1FA27@G11@XFT28#ZPFT27ZP1024FP610D@524D!FO30@SFL8FE22@ |
|||
[===================================================================] |
|||
[The following, without the comments and white space, might have |
|||
been input from a separate tape.] |
|||
E25K TM GK |
|||
E19Z [define entry point] |
|||
PF [acc = 0 on entry] |
|||
[Integers supplied by user: (1) seed for LCG; (2) list of numbers of trials |
|||
for which to print result; increasing order, terminated by 0. |
|||
To be read by library subroutine R4; sign comes after value.] |
|||
987654321+100+1000+10000+100000+0+ |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
TRIALS EST PI/10 |
|||
100 0.32400 |
|||
1000 0.31319 |
|||
10000 0.31371 |
|||
100000 0.31410 |
|||
</pre> |
|||
=={{header|Elixir}}== |
=={{header|Elixir}}== |