理想及其算法設(shè)計

宗 慧，曹子寧，王?？?/p>

（南京航空航天大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院，江蘇南京 211106）

隨機(jī)數(shù)在統(tǒng)計理論、科學(xué)研究和工程設(shè)計中具有重要作用［1］。盡管在真實的物理世界中可以得到隨機(jī)數(shù)，但限于生成速度和生成方法，真隨機(jī)數(shù)的獲取極為困難。隨著電子數(shù)字計算機(jī)的問世，采用算法產(chǎn)生隨機(jī)數(shù)成為一種行之有效的方法，所產(chǎn)生的這些數(shù)字稱作偽隨機(jī)數(shù)。當(dāng)今隨機(jī)數(shù)在現(xiàn)代科學(xué)中的應(yīng)用非常廣泛，如計算科學(xué)、通信、網(wǎng)絡(luò)安全、雷達(dá)探測、統(tǒng)計分析、數(shù)值模擬等領(lǐng)域都有大量的應(yīng)用［2－3］。 然而，偽隨機(jī)數(shù)發(fā)生器的質(zhì)量對這些應(yīng)用的成效至關(guān)重要。因此，偽隨機(jī)數(shù)發(fā)生器質(zhì)量受到科技界的關(guān)注。德國Federal Office for Information Security （Bundesamt für Sicherheit in der Informationstechnik，BSI）評估隨機(jī)生成器質(zhì)量的主要準(zhǔn)則要求所產(chǎn)生的序列的相同概率低，符合統(tǒng)計學(xué)的平均性。美國國家標(biāo)準(zhǔn)與技術(shù)研究院（National Iustitute of Standardsand Technology，NIST）側(cè)重于數(shù)學(xué)性能而制定了隨機(jī)性測試方法，主要包括對頻數(shù)、塊內(nèi)頻數(shù)、矩陣秩、離散傅里葉變換以及幾種組合的測試［4］。 根據(jù) BSI主要準(zhǔn)則，文獻(xiàn)［5］利用 Monte Carlo計算值的方法建立了一種評價偽隨機(jī)數(shù)發(fā)生器質(zhì)量的方法，該方法將理想π值作為標(biāo)準(zhǔn)，以實際計算出的值與其之差評價了多種偽隨機(jī)數(shù)發(fā)生器的質(zhì)量。

1 面積之比——計算的基本思想

因此

含有內(nèi)切圓的正方形如圖1（a）所示。如果該正方形由N2個面積為1的小正方形組成，那么正方形的面積就等于所有小正方形的面積之和，即As＝N2。

如果把圖1（a）中的每個小正方形用點來替代，如圖1（b）所示，那么點的數(shù)量之和就等于正方形的面積，如果把點看作單位1，當(dāng)點足夠小，且當(dāng)相鄰的兩點的距離趨于零時，那么區(qū)域內(nèi)點的數(shù)量之和就可以等價于該區(qū)域的面積，以此類推，根據(jù)式（1）知：計算值時，只需要計算內(nèi)切圓面積Ac與正方形面積As之比。因此，可以先將圖1（a）中的所有小正方形縮為圖1（b）中的點，當(dāng)相鄰兩點的距離趨于零時，就可以用內(nèi)切圓和正方形內(nèi)的點分別代替式（1）中的Ac和As進(jìn)行計算。這就可以在程序中方便地以點代面實現(xiàn)對值的計算，從而避免了程序從面的角度計算π值的困境。

圖1 以點代面的原理

定義1給定一個由均勻分布的充分多的點組成的含有內(nèi)切圓的正方形，并設(shè)正方形內(nèi)的點數(shù)為Xs，內(nèi)切圓內(nèi)的點數(shù)為 Xc，有

定理1

為了使證明簡潔，先回憶一下數(shù)學(xué)中對點、線和數(shù)軸的有關(guān)描述。

點是沒有部分的，即它只表示位置而無大小。兩個不同的點可以確定一條直線［15］。如果在直線上確定一個點O為原點，并指定一個方向為正向，則該直線為數(shù)軸，顯然數(shù)軸是由點的集合構(gòu)成的。數(shù)軸上的每個點當(dāng)且僅當(dāng)對應(yīng)一個實數(shù)。由于實數(shù)是稠密的，即任意兩個不同的實數(shù)間必存在另一個實數(shù)，所以數(shù)軸上的點是稠密的。又由于實數(shù)集是無限集，所以實數(shù)上的點集也是無限集。

證明：根據(jù)定義1，設(shè)含有內(nèi)切圓的正方形的邊長為R，并構(gòu)成如圖2所示的直角坐標(biāo)。此時的3個點（0，0），（0，R）和（R，0）就確定了一個平面。 在向該平面內(nèi)的正方形投放均勻分布的充分多的點時，就可以由點數(shù)計算值。然而，只有在這些投放的點是稠密的，即是無限時，才能得到與理論一致的值，即

圖2 邊長為R的正方形所在的平面

利用式（2）進(jìn)行計算，正方形內(nèi)必須具有一定的點數(shù)，才能得到滿足所要求精度的′。例如，正方形內(nèi)只有100個點，顯然，此時計算出來的′與理論值的誤差很大，因此使用較少的點取代面積來計算′沒有意義。那么在正方形內(nèi)投放多少個點才是有意義的呢？這是本文要解決的關(guān)鍵問題。下面就此進(jìn)行討論。

定義2給定一個由102r個均勻分布的點組成的含有內(nèi)切圓的正方形。稱圓內(nèi)的點數(shù)為理想點數(shù)，記作XI；由該圖形求出的′稱作理想，記作I，且

定理2給定一個由102r個均勻分布的點組成的正方形，并設(shè)I與π的計算誤差為10－m，其中m＞0，有

證明：（i）假設(shè)落入正方形內(nèi)點數(shù)為102r，并設(shè)

即

根據(jù)式（4），有

即

兩邊取對數(shù)，有

故

根據(jù)式（4），有

即

兩邊取對數(shù)，可得

故

根據(jù)（a）和（b），式（5）得證。

故m ＜0，與定理2所給條件矛盾。

（iii）的證明與（i）類似，略。

定理2表明了r與m之間的關(guān)系。根據(jù)定理2可知，當(dāng)正方形內(nèi)點的數(shù)量給定，那么I與的計算誤差就可以得到，同時I的精度也就確定了。

3 計算πI的算法設(shè)計

為了減少計算量，在算法設(shè)計時取如圖3（a）所示的含四分之一內(nèi)切圓的正方形計算I值，首先在正方形內(nèi)均勻地投放n2個點，這里的 n＝10r；然后統(tǒng)計內(nèi)切圓內(nèi)點的總數(shù)，最后根據(jù)式（4）計算I。

（count＝0） ∥count為內(nèi)切圓內(nèi)點數(shù)之和

在這個算法中采取的是利用循環(huán)結(jié)構(gòu)依次統(tǒng)計內(nèi)切圓內(nèi)的所有點，如果投放的點數(shù)非常多，那么統(tǒng)計的點的數(shù)量也會加大，從而導(dǎo)致運(yùn)算的時間變長，因此本文對這個算法進(jìn)行了改進(jìn)以縮短計算時間。如圖3（a）所示，可以統(tǒng)計正方形內(nèi)內(nèi)切圓外的灰色點的數(shù)量，記為Yc，這個點的數(shù)量是明顯少于內(nèi)切圓內(nèi)的點的數(shù)量Xc的，所以同樣投點的情況下，計算的次數(shù)減少了，計算的時長也縮短了，但是這種算法減少的計算次數(shù)還不是特別顯著。在此基礎(chǔ)上本文又做了一個改進(jìn)，依然是統(tǒng)計 Yc，如圖 3（b）所示，正方形的對角線與內(nèi)切圓的交點位于坐標(biāo)（0.707 106 78n， 0.707 106 78n）處，因此只需要統(tǒng)計S1區(qū)域內(nèi)的點數(shù) S1＿count，而 S2區(qū)域內(nèi)的點數(shù)是可以根據(jù)直接計算得到，記為S2＿count，然后可以得到 Yc＝ 2S1＿count＋S2＿count，按照這個思想設(shè)計的算法運(yùn)算的次數(shù)大幅度降低，計算的時長也縮短了。

圖3 計算πI的算法設(shè)計

具體算法如下：

（S1＿count＝0） ∥S1＿count為內(nèi)切圓外S1區(qū)域內(nèi)點數(shù)之和

4 實驗結(jié)果與分析

依據(jù)表1可以看出：

表 1 不同的投點數(shù)對應(yīng)的I值（π＝3.141 592 653 6……）

表 1 不同的投點數(shù)對應(yīng)的I值（π＝3.141 592 653 6……）

注：列A表示IDEAL＿1算法實驗結(jié)果，列B表示IDEAL＿2算法實驗結(jié)果

總投點數(shù) r m πI 計算時間／ms A B A B A B 25 000 000 3.698 970 004 3 3.093 9 3.470 8 3.140 787 040 0 3.141 254 400 0 18 4 100 000 000 4.000 000 000 0 3.395 6 3.774 3 3.141 190 520 0 3.141 424 520 0 72 15 2 500 000 000 4.698 970 004 3 4.095 6 4.095 5 3.141 512 417 6 3.141 512 401 6 1 702 394 10 000 000 000 5.000 000 000 0 4.396 8 4.396 7 3.141 552 545 6 3.141 552 541 6 6 795 1 561 250 000 000 000 5.698 970 004 3 5.095 9 5.095 9 3.141 584 635 8 3.141 584 635 4 169 707 39 047 1 000 000 000 000 6.000 000 000 0 5.397 5 5.397 5 3.141 588 649 6 3.141 588 649 3 679 234 156 640 25 000 000 000 000 6.698 970 004 3 6.096 8 6.478 9 3.141 591 853 3 3.141 592 321 6 16 794 482 3 858 138 100 000 000 000 000 7.000 000 000 0 6.397 8 6.779 7 3.141 592 253 5 3.141 592 487 5 67 405 604 15 442 081

表2 不同有效位的理想值對應(yīng)的最低投點數(shù)

表2 不同有效位的理想值對應(yīng)的最低投點數(shù)

總投點數(shù) r 圓內(nèi)點數(shù) πI m 有效位6 512 704 3.406 9 5 112 486 3.140 008 205 5 2.800 1 2 45 832 900 3.830 6 35 990 287 3.141 000 198 5 3.227 3 3 1 873 158 400 4.636 3 1 471 131 781 3.141 500 005 6 4.033 2 4 2 274 380 691 025 6.178 4 1 471 131 782 3.141 590 000 0 5.576 2 5 37 471 488 988 816 6.786 9 29 430 051 739 428 3.141 592 000 0 6.184 7 6

最后本文結(jié)合表1和表2中的數(shù)據(jù)，繪制了r和m之間的關(guān)系，如圖4所示。根據(jù)該圖可以得到m值隨著r值的增長而增長。圖4中的曲線顯示，對于投點數(shù)（10r）2和精度（10－m）之間的關(guān)系，IDEAL＿1和 IDEAL＿2兩種算法的結(jié)果是一致的。

圖4 投點數(shù)（10r）2和精度（10－m）之間的關(guān)系

5 結(jié)束語

也可以用其他需要用到隨機(jī)數(shù)發(fā)生器計算的常數(shù)或諸如面積、體積之類的方法去評價偽隨機(jī)數(shù)發(fā)生器，比如e，這是一種比較簡便的方法。