Monte Carlo methods: Difference between revisions

← Older edit

Monte Carlo methods (view source)

Revision as of 10:48, 4 January 2024

6,469 bytes added , 4 months ago

m

→‎{{header|Wren}}: Minor tidy

PureFox

9,476

edits

Revision as of 05:01, 1 May 2023 (view source) Elfish (talk \| contribs) No edit summary ← Older edit		Latest revision as of 10:48, 4 January 2024 (view source) PureFox (talk \| contribs) m (→‎{{header\|Wren}}: Minor tidy)
(8 intermediate revisions by 5 users not shown)
Line 1,022: =={{header\|EasyLang}}== <syntaxhighlight lang="text">~~func mc n . .~~ func mc n . ~~for i range n~~ for xi = ~~randomf~~1 to n y x = randomf if x * x + y * y <= 1randomf ~~hit~~if x * x += y * y < 1 . hit += 1 . . ~~print 4.0 * hit / n~~ return 4.0 * hit / n . numfmt 4 0 ~~fscale 4~~ ~~call~~print mc 10000 ~~call~~print mc 100000 ~~call~~print mc 1000000 ~~call~~print mc 10000000~~</syntaxhighlight>~~ </syntaxhighlight> Output: 3.1292 Line 1,042 ⟶ 1,044: 3.1407 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}}== Line 1,603 ⟶ 1,733: piMC 10^i.7 4 2.8 3.24 3.168 3.1432 3.14256 3.14014</syntaxhighlight> '''Alternative Tacit Solution:''' <syntaxhighlight lang="j">pimct=. (4 +/ % #)@:(1 >: \|)@:(? j. ?)@:($&0)"0 (,. pimct) 10 ^ 3 + i.6 1000 3.168 10000 3.122 100000 3.13596 1e6 3.1428 1e7 3.14158 1e8 3.14154</syntaxhighlight> =={{header\|Java}}== Line 1,865 ⟶ 2,006: for n in 10 .^ (3:8) p = monteπ(n) println("$(lpad(n, 9)): π ≈ $(lpad(p, 10)), pct.err = ", @sprintf("%2.5f%%", 100 * abs(p - π) / π)) end</syntaxhighlight> {{out}} <pre> 1000: π ≈ 3.224, pct.err = 0.02623% 10000: π ≈ 3.1336, pct.err = 0.~~00254~~254% 100000: π ≈ 3.13468, pct.err = 0.~~00220~~220% 1000000: π ≈ 3.14156, pct.err = 0.~~00001~~001% 10000000: π ≈ 3.1412348, pct.err = 0.~~00011~~011% 100000000: π ≈ 3.14123216, pct.err = 0.~~00011~~011%</pre> =={{header\|K}}== Line 2,688 ⟶ 2,829: monteCarlo(10000) = 3.00320000 monteCarlo(100000) = 2.99536000 </pre> =={{header\|RPL}}== ≪ 0 1 3 PICK '''START''' RAND SQ RAND SQ + 1 < + '''NEXT''' SWAP / 4 * ≫ '<span style="color:blue">MCARL</span>' STO 100 <span style="color:blue">MCARL</span> 1000 <span style="color:blue">MCARL</span> 10000 <span style="color:blue">MCARL</span> 100000 <span style="color:blue">MCARL</span> {{out}} <pre> 4: 3.2 3: 3.084 2: 3.1684 1: 3.14154 </pre> Line 3,058 ⟶ 3,219: {{trans\|Kotlin}} {{libheader\|Wren-fmt}} <syntaxhighlight lang="~~ecmascript~~wren">import "random" for Random import "./fmt" for Fmt var rand = Random.new()