Monte Carlo methods: Difference between revisions
m
→{{header|Wren}}: Minor tidy
(Monte Carlo methods in various BASIC dialents (BASIC256, Run BASIC, True BASIC and Yabasic)) |
m (→{{header|Wren}}: Minor tidy) |
||
| (10 intermediate revisions by 7 users not shown) | |||
Line 931:
The script give in reality an output formatted in HTML
<pre>π ≅ 3.14139</pre>
=={{header|Delphi}}==
{{works with|Delphi|6.0}}
{{libheader|SysUtils,StdCtrls}}
<syntaxhighlight lang="Delphi">
function MonteCarloPi(N: cardinal): double;
{Approximate Pi by seeing if points fall inside circle}
var I,InsideCnt: integer;
var X,Y: double;
begin
InsideCnt:=0;
for I:=1 to N do
begin
{Random X,Y = 0..1}
X:=Random;
Y:=Random;
{See if it falls in Unit Circle}
if X*X + Y*Y <= 1 then Inc(InsideCnt);
end;
{Because X and Y are squared, they only fall with 1/4 of the circle}
Result:=4 * InsideCnt / N;
end;
procedure ShowOneSimulation(Memo: TMemo; N: cardinal);
var MyPi: double;
begin
MyPi:=MonteCarloPi(N);
Memo.Lines.Add(Format('Samples: %15.0n Pi= %2.15f',[N+0.0,MyPi]));
end;
procedure ShowMonteCarloPi(Memo: TMemo);
begin
ShowOneSimulation(Memo,1000);
ShowOneSimulation(Memo,10000);
ShowOneSimulation(Memo,100000);
ShowOneSimulation(Memo,1000000);
ShowOneSimulation(Memo,10000000);
ShowOneSimulation(Memo,100000000);
end;
</syntaxhighlight>
{{out}}
<pre>
Samples: 1,000 Pi= 3.156000000000000
Samples: 10,000 Pi= 3.152000000000000
Samples: 100,000 Pi= 3.142920000000000
Samples: 1,000,000 Pi= 3.140864000000000
Samples: 10,000,000 Pi= 3.141990800000000
Samples: 100,000,000 Pi= 3.141426720000000
</pre>
=={{header|E}}==
Line 965 ⟶ 1,022:
=={{header|EasyLang}}==
<syntaxhighlight lang="text">
func mc n .
for
.
.
return 4.0 * hit / n
.
numfmt 4 0
</syntaxhighlight>
Output:
3.1292
Line 985 ⟶ 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,546 ⟶ 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,808 ⟶ 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.
100000: π ≈ 3.13468, pct.err = 0.
1000000: π ≈ 3.14156, pct.err = 0.
10000000: π ≈ 3.1412348, pct.err = 0.
100000000: π ≈ 3.14123216, pct.err = 0.
=={{header|K}}==
Line 2,631 ⟶ 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 2,840 ⟶ 3,058:
1000000: 3.142344
</pre>
=={{header|SparForte}}==
As a structured script.
<syntaxhighlight lang="ada">#!/usr/local/bin/spar
pragma annotate( summary, "monte" )
@( description, "A Monte Carlo Simulation is a way of approximating the" )
@( description, "value of a function where calculating the actual value is" )
@( description, "difficult or impossible. It uses random sampling to define" )
@( description, "constraints on the value and then makes a sort of 'best" )
@( description, "guess.'" )
@( description, "" )
@( description, "Write a function to run a simulation like this with a" )
@( description, "variable number of random points to select. Also, show the" )
@( description, "results of a few different sample sizes. For software" )
@( description, "where the number pi is not built-in, we give pi to a couple" )
@( description, "of digits: 3.141592653589793238462643383280 " )
@( see_also, "http://rosettacode.org/wiki/Monte_Carlo_methods" )
@( author, "Ken O. Burtch" );
pragma license( unrestricted );
pragma restriction( no_external_commands );
procedure monte is
function pi_estimate (throws : positive) return float is
inside : natural := 0;
begin
for throw in 1..throws loop
if numerics.random ** 2 + numerics.random ** 2 <= 1.0 then
inside := @ + 1;
end if;
end loop;
return 4.0 * float (inside) / float (throws);
end pi_estimate;
begin
? " 1_000:" & strings.image (pi_estimate ( 1_000))
@ " 10_000:" & strings.image (pi_estimate ( 10_000))
@ " 100_000:" & strings.image (pi_estimate ( 100_000))
@ " 1_000_000:" & strings.image (pi_estimate ( 1_000_000));
end monte;</syntaxhighlight>
=={{header|Stata}}==
Line 2,960 ⟶ 3,219:
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
import "./fmt" for Fmt
var rand = Random.new()
| |||