我在 PSTricks 中编写了一个宏，利用所谓的“蒙特卡罗方法”可以近似地计算出圆周率 (π)。
将设置一个表格，其中包含两个随机选择的数字 x 和 y（每个数字介于 0 和 1 之间）。它们每个都有一个单独的计数器，取决于 x^2 + y^2 > 1 还是 x^2 + y^2 <= 1 现在我想在 LaTeX/L3 中重写它。有人可以给我提示如何开始吗？非常感谢。
\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}
答案1
例如：
\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}
答案的细化：
\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
仅使用 TeX 整数就可以获得更快的执行时间。
\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}
但是仅使用随机整数有一个警告<32768……：
代码：
\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}
然后，我想到使用随机整数<2**28而不是<2**15。请注意，PDFTeX RNG 内部使用 28 位整数。关键是下面的代码仍然只使用 TeX 整数，但它利用了\numexpr临时切换到双精度的能力。
因此，该代码并不比之前的代码慢很多。
\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}
得出的结果为：
我不知道它（非常缓慢地）收敛到的极限的确切值，但它会比3.141713824123144 . . .上面的更接近 Pi。
我运行了“大”代码，并使用了一个输出到文件的变体。在试验结束时10**9（等待时观看了《飞出个未来》一集），我得到了这个
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
与我在相同情况下使用随机整数<32768得到的“小”代码相比，使用随机整数的结果要好得多。3.14170127610**93.141574576< 268435456
看起来“大”代码的执行只是50%比“小”代码慢一些。
答案3
这是另一种方法，可以通过 JavaScript 在 PDF 文件中以交互方式完成计算。但是，需要 Adobe Reader 才能查看文档：
\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}
