Rewrite a postscript macro in LaTeX
E. g.:
\documentclass{article}
\usepackage{expl3}
\ExplSyntaxOn
\newcommand*\piMC[1]{
\pi_mc:n{#1}
}
\cs_new_nopar:Nn\pi_mc:n{
\int_zero_new:N\l_M_int
\prg_replicate:nn{#1}{
\fp_compare:nNnT{rand()**2+rand()**2}<{1}{\int_incr:N\l_M_int}
}
\fp_eval:n{4*\int_use:N\l_M_int/#1}
}
\ExplSyntaxOff
\begin{document}
\begin{tabular}{rl}
$N$ & $\pi$ \\\hline
10 & \piMC{10}\\
100 & \piMC{100}\\
1000 & \piMC{1000}\\
10000 & \piMC{10000}
\end{tabular}
\end{document}
A refinement of the answer:
\documentclass{article}
\usepackage{xparse}
\ExplSyntaxOn
\NewDocumentCommand\mcpi {m}
{
\thomas_mcpi:n { #1 }
}
\int_new:N \l_thomas_mcpi_in_int
\cs_new_protected:Nn \thomas_mcpi:n
{
\int_zero:N \l_thomas_mcpi_in_int
\prg_replicate:nn { #1 }
{
\fp_compare:nNnT { rand()**2+rand()**2 }<{1}
{ \int_incr:N \l_thomas_mcpi_in_int }
}
\fp_eval:n { 4*\l_thomas_mcpi_in_int/#1 }
}
\ExplSyntaxOff
\begin{document}
\begin{tabular}{rl}
$N$ & $\pi$ \\\hline
10 & \mcpi{10}\\
100 & \mcpi{100}\\
1000 & \mcpi{1000}\\
10000 & \mcpi{10000}
\end{tabular}
\end{document}
Much much faster execution times can be obtained from using only TeX integers.
\documentclass{article}
\usepackage{xintkernel}% for xintreplicate (which is a clone of the expl3 replicate)
\newcount\montecarlocount
\newcount\tmpcnta
\newcount\tmpcntb
\newcount\twotothefifteen
\twotothefifteen "8000 % 8*16**3 = 2**15 = 32768
\newcount\twotothethirty
\twotothethirty "40000000 % 4*16**7 = 2**30 = 1073741824
\makeatletter
\newcommand\onestep{%
\tmpcnta\pdfuniformdeviate\twotothefifteen\relax
\tmpcntb\pdfuniformdeviate\twotothefifteen\relax
% A**2 + B**2 < R**2 ?
\ifnum\numexpr\tmpcnta*\tmpcnta+\tmpcntb*\tmpcntb < \twotothethirty
\advance\montecarlocount \@ne
\fi
}%
\newcommand\piMC[1]{%
\montecarlocount \z@
\romannumeral\xintreplicate{#1}\onestep
\strip@pt\dimexpr\numexpr4*\montecarlocount*65536/#1 sp\relax
}
\newcount\gtotalcount
\newcount\gmontecarlocount
\newcommand\piMCstart{%
\global\gtotalcount\z@
\global\gmontecarlocount\z@
}
\newcommand\piMCmore[2]{%
\montecarlocount \z@
\romannumeral\xintreplicate{#1}{\romannumeral\xintreplicate{#2}\onestep}%
\global\advance\gtotalcount \numexpr#1*#2\relax
\global\advance\gmontecarlocount\montecarlocount
\strip@pt\dimexpr\numexpr4*\gmontecarlocount*65536/\gtotalcount sp\relax
}
\makeatother
\begin{document}
\pdfsetrandomseed 12345 % nota bene: automatically global
\begin{tabular}{rlrr}
$N$ & $\pi$ & hits & total\\\hline
10 & \piMC{10}\\
100 & \piMC{100}\\
1000 & \piMC{1000}\\
10000 & \piMC{10000}\\
100000 & \piMCstart
\piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
% \romannumeral\xintreplicate{8}{+100000 & \piMCmore{100}{1000}
% &\the\gmontecarlocount&\the\gtotalcount\\}%
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount
\end{tabular}
\end{document}
There is a caveat though with using only random integers <32768...:
Code:
\documentclass{article}
\usepackage{amssymb}
\usepackage{xintexpr}
\newcount\indicatorcount
\begin{document}
% we count exactly all pairs (x,y) with
% 1) 0<= x < 32768, 0 <= y< 32768
% 2) x**2+y**2 < 2**30 = 1073741824
% for this given x we need to count how many
% y < 32768 are there with y**2 < Z = 2**30 - x**2
% This means that y**2 is at most Z - 1
% This means that y is at most sqrt(Z-1)
% This means that y is at most trunc(sqrt(Z-1)) (automatically < 32768)
% This means that there are exactly 1 + trunc(sqrt(Z-1)) solutions
% for each x.
% we can use sqrt() function in \xintiiexpr to get
% the truncated square root.
% We use \xintexpr syntax to do the computation in one line
\edef\zzz{\xintiiexpr 2**30-1\relax}% don't use \xinttheiiexpr, keep it in
% private, so faster
\edef\www{\xintiiexpr 2**15\relax}
\edef\totalnumberofsolutions
{\xinttheiiexpr iter(0; (i=\www)?{abort}
{@ + 1 + sqrt(\zzz - i**2)}
, i = 0++)\relax}
\totalnumberofsolutions\space pairs $(x, y)$ exist with verify
\begin{enumerate}
\item $0\leqslant x, y < 32768$,
\item $x^2 +y^2 < 32768^2$.
\end{enumerate}
This means that Monte Carlo methods using only \verb|\pdfuniformdeviate 32768|
will (slowly) give an approximation, not of $\pi$ but of:
% 4*\totalnumberofsolutions would go in \numexpr and create arithmetic overflow!
% so we can not do this: \xintTrunc{15}{4*\totalnumberofsolutions/1073741824}
% 268435456 = 1073741824/4
% so we do that:
\[\xintTrunc{15}{\totalnumberofsolutions/268435456}\dots\]
\end{document}
Then, I figured out to use random integers <2**28 rather than <2**15. Note that internally the PDFTeX RNG works with 28bits integers. The point is that the code below still uses only TeX integers, but it exploits \numexpr ability to temporary switch to doubled-precision.
Thus, the code is not much slower than the earlier one.
\documentclass{article}
\usepackage{xintkernel}% for xintreplicate (which is a clone of the expl3 replicate)
\newcount\montecarlocount
\newcount\tmpcnta
\newcount\tmpcntb
\newcount\tmpcntc
\newcount\tmpcntd
\newcount\maxrand
\newcount\threemaxrand
%% 16**7 = 2**28. The pdftex RNG works natively with 28bits integers
\maxrand "10000000
\threemaxrand "30000000
% I considered using \maxint for increased precision in the comparison below
% but with \maxrand gives same result
\newcount\maxint
% % % 2**31 - 1 (8 times \cntc minus 1)
\maxint "7FFFFFFF
\makeatletter
\newcommand\onestep{%
\tmpcnta\pdfuniformdeviate\maxrand\relax % A
\tmpcntb\pdfuniformdeviate\maxrand\relax % B
% \ifnum\tmpcnta<\tmpcntb % this would make A >= B in what follows
% \tmpcntc\tmpcntb
% \tmpcntb\tmpcnta
% \tmpcnta\tmpcntc
% \fi % but we don't use it finally
% This is commented-out too, it compares A+B with 1.5M
% \tmpcntc\tmpcnta
% \advance\tmpcntc\tmpcntb % A + B
% \multiply\tmpcntc \tw@
% \ifnum\tmpcntc<\threemaxrand
\tmpcntd\maxrand % M
\advance\tmpcntd-\tmpcntb % D = M - B
\ifnum\tmpcntd>\tmpcnta %% always true if A = 0 !!
\advance\montecarlocount \@ne % M - B > A => A^2 + B^2 < M^2
\else
% here M <= A + B
\tmpcntc\tmpcntb
\advance\tmpcntc\maxrand % C = M + B
\ifnum % use temporary double precision!
% I considered using \maxint for increased precision but does not change result
\numexpr\maxrand*\tmpcnta/\tmpcntc % A/(M + B) (is < 1)
<\numexpr\maxrand*\tmpcntd/\tmpcnta\relax % (M - B)/A (is <= 1)
% the comparison will not be very precise when both A<<M+B and M-B<<A
% this can happen only if M-B<<M+B, i.e. for B close to M.
%
% If B is arranged to be <= A, both A and B are then close to M and A+B is close to
% 2M. To make sure this does not happen we could require A+B < 1.5 M
% (if A+B>=1.5M, then A^2+B^2>M^2 ...)
%
% But in fact even without precautions same result was obtained...
% .. and we can gain 20% time be suppressing all precautions...
% seems that comparison above with 30bits is precise enough
% So I have commented out the code which made sure B<= A and A+B< 1.5M
\advance\montecarlocount \@ne
\fi
\fi
% \fi
}%
\newcommand\piMC[1]{%
\montecarlocount \z@
\romannumeral\xintreplicate{#1}\onestep
\strip@pt\dimexpr\numexpr4*\montecarlocount*65536/#1 sp\relax
}
\newcount\gtotalcount
\newcount\gmontecarlocount
\newcommand\piMCstart{%
\global\gtotalcount\z@
\global\gmontecarlocount\z@
}
\newcommand\piMCmore[2]{%
\montecarlocount \z@
\romannumeral\xintreplicate{#1}{\romannumeral\xintreplicate{#2}\onestep}%
\global\advance\gtotalcount \numexpr#1*#2\relax
\global\advance\gmontecarlocount\montecarlocount
\strip@pt\dimexpr\numexpr4*\gmontecarlocount*65536/\gtotalcount sp\relax
}
\makeatother
\begin{document}
\pdfsetrandomseed 12345 % nota bene: automatically global
\begin{tabular}{rlrr}
$N$ & $\pi$ & hits & total\\\hline
10 & \piMC{10}\\
100 & \piMC{100}\\
1000 & \piMC{1000}\\
10000 & \piMC{10000}\\
100000 & \piMCstart
\piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
% \romannumeral\xintreplicate{38}{+100000 & \piMCmore{100}{1000}
% &\the\gmontecarlocount&\the\gtotalcount\\}%
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount\\
+100000 & \piMCmore{100}{1000}&\the\gmontecarlocount&\the\gtotalcount
\end{tabular}
\end{document}
This produces:
I don't know the exact value of the limit to which it will (very slowly) converges, but it will be closer to Pi than the 3.141713824123144 . . . above.
I ran the "big" code, with a variant which outputs to a file. At the end of the 10**9 trials (watch a Futurama episode while waiting) I got this
999990000: 3.141574507745077
999991000: 3.141574510170592
999992000: 3.141574484595877
999993000: 3.141574483021381
999994000: 3.141574581447489
999995000: 3.141574591872959
999996000: 3.141574582298329
999997000: 3.141574628723886
999998000: 3.141574595149190
999999000: 3.141574585574586
1000000000: 3.141574576000000
In contrast with the "small" code (using random ints <32768) I had gotten under same circumstances 3.141701276 after 10**9 random points. Clearly the 3.141574576 is much better from using random ints < 268435456.
It seems execution of the "big" code is only about 50% slower than "small" code.
Here another approach to get the calculations interactively done by JavaScript within the PDF file. However one needs Adobe Reader to view the document:
\documentclass{article}
\usepackage{eforms}
\begin{defineJS}{\myPi}
var r = 5;
var points_total = 0;
var points_inside = 0;
var rowA ="";
var f1 = this.getField("PiField");
while(1){
points_total++;
var x = Math.random()*r*2-r;
var y = Math.random()*r*2-r;
if (Math.pow(x,2) + Math.pow(y,2) < Math.pow(r,2))
points_inside++;
if (points_total \% 100 == 0)
util.printf("\%,2 .9f",(4*points_inside/points_total))+"\\r";
if (points_total == 1000)
break;
var row="Total "
+util.printf("\%,2 d",points_total)
+": pi = "
+util.printf("\%,2 .9f",(4*points_inside/points_total))+"\\r";
rowA+=row;
}
f1.value=rowA;
\end{defineJS}
\begin{document}
\pushButton[\CA{Calculate Pi}\A{\JS{\myPi}}]{MCpi}{}{11bp} \qquad \pushButton[\CA{Reset}\A{\JS{this.resetForm("PiField");}}]{reset}{1.5cm}{11bp}
\bigskip
\textField[\Ff{\FfMultiline}]{PiField}{8cm}{11bp*40}
\end{document}
