8

I wrote a macro in PSTricks with which the number pi is calculated approximately with the so-called "Monte-Carlo-Method".

A table will be setup with two randomly chosen numbers x and y (each of them between 0 and 1). Each of them has a seperate counter depending if x^2 + y^2 > 1 or x^2 + y^2 <= 1 Now I'd like to rewrite that in LaTeX/L3. Is there anybody to give me a hint how to start? Many thanks.

\documentclass[a4paper,dvipsnames,svgnames]{article}
\usepackage[margin=1.0cm]{geometry}
\usepackage{pstricks,pstricks-add}

\makeatletter
\define@key[psset]{}{colorF}[blue]{\pst@getcolor{#1}\psk@colorF}
\define@key[psset]{}{colorT}[red]{\pst@getcolor{#1}\psk@colorT}
\define@key[psset]{}{colorPi}[green]{\pst@getcolor{#1}\psk@colorPi}
\psset{colorF=blue,colorT=red,colorPi=green}
%
\def\psRandomPiTable{\pst@object{psRandomPiTable}}
\def\psRandomPiTable@i#1{%
 \pst@killglue
  \addbefore@par{colorF=DarkBlue,colorT=BrickRed,colorPi=Green,fontscale=7}%
  \begin@SpecialObj
  \addto@pscode{
        /s1 { /Symbol findfont \psk@fontscale\space 4 add scalefont setfont } bind def
        /s2 { \psk@PSfont findfont \psk@fontscale\space scalefont setfont } bind def
        /s3 { \psk@PSfont findfont \psk@fontscale\space 4 add scalefont setfont } bind def
        /zaehlerT 0 def
        /zaehlerF 0 def
        /dec 3 def
        /ru \psk@fontscale\space 0.5 add def
        /re 30 def
        realtime srand Rand pop
    4 dict begin
       1 1 #1 {
                /i ED
                /x Rand def
                /y Rand def
                /r { x dup mul y dup mul add sqrt } bind def
                s2
                /showvalue { 10 dec exp mul round 10 dec exp div 9 string cvs dot2comma show } bind def
                0 i 1 add ru mul neg moveto %
                i 3 string cvs show
                re i 1 add ru mul neg moveto %
                x showvalue
                2 re mul i 1 add ru mul neg moveto %
                y showvalue
                gsave
                r 1 lt
                    { /zaehlerF zaehlerF 1 add def \pst@usecolor\psk@colorF 3 re mul i 1 add ru mul neg moveto }
                    { /zaehlerT zaehlerT 1 add def \pst@usecolor\psk@colorT 4 re mul i 1 add ru mul neg moveto }
                ifelse %
                r showvalue
                r 1 lt
                    { 5 re mul i 1 add ru mul neg moveto (Nr. ) show zaehlerF 10 string cvs show }
                    { 5 re mul i 1 add ru mul neg moveto (Nr. ) show zaehlerT 10 string cvs show }
                ifelse %
                grestore
               } for
        gsave \pst@usecolor\psk@colorF \psk@PSfont findfont \psk@fontscale\space 4 add scalefont setfont
        2 re mul 0.4 ru mul neg moveto (Summe: ) show zaehlerF 10 string cvs show
        grestore
        gsave \pst@usecolor\psk@colorT\psk@PSfont findfont \psk@fontscale\space 4 add scalefont setfont
        4 re mul 0.4 ru mul neg moveto (Summe: ) show zaehlerT 10 string cvs show
        grestore
        gsave \pst@usecolor\psk@colorPi
        0 re mul 0.4 ru mul neg moveto s1 (\string\160 \string\273 ) show s3 4 zaehlerF i div mul 10 string cvs dot2comma show
        grestore
    end
  }%
  \end@SpecialObj
 \ignorespaces
}

\begin{document}
\section*{$\pi$-Monte-Carlo}
\psRandomPiTable{100}
\end{document}

enter image description here

  • It would be easier to reimplement this if you provided to original (infix) formulae. – AlexG May 3 '18 at 11:53
  • 1
    I think I got the idea, by reading here. – AlexG May 3 '18 at 11:58
  • @AlexG yes, that is it, in the code: 4 zaehlerF i div mul means 4*Z/N, in the example: 4*79/(79+21)=3,16 – user139826 May 3 '18 at 12:02
10

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}
  • 2
    I took the liberty of adding a polished version of your code, more adherent to the guidelines. Note that \thomas_mcpi:n (your \pi_mc:n) should be “protected”, because it does an assignment; the int variable is better allocated outside the function. – egreg May 3 '18 at 16:43
  • @egreg Within a package, should any macro doing assignments be defined \protected, even those defined for internal use, or just those exposed to users of the package? – AlexG May 4 '18 at 8:44
  • @AlexG Yes, it should – egreg May 4 '18 at 9:29
  • Ok, thank you, @egreg. But then, files like latex.ltx should be full of \protected. – AlexG May 4 '18 at 9:31
  • 1
    @AlexG \protected wasn't normal back when latex.ltx was written. If it was written nowadays it would probably have many many \protected inside. – Manuel May 4 '18 at 12:28
5

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}

enter image description here

There is a caveat though with using only random integers <32768...:

enter image description here

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:

enter image description here

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.

  • you once again thought a lot and were very successful. Thank you very much. – user139826 May 4 '18 at 11:47
5

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}

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy