基于FPGA的高斯分布隨機(jī)數(shù)的生成*

姚若河 徐新才

(華南理工大學(xué) 電子與信息學(xué)院，廣東 廣州 510640)

蒙特卡羅模擬廣泛應(yīng)用于數(shù)學(xué)、計算物理學(xué)、工程技術(shù)以及經(jīng)濟(jì)金融等方面［1-5］，用以計算多重積分、解積分方程和隨機(jī)微分方程等;或者直接模擬隨機(jī)過程(布朗運(yùn)動、動物的生態(tài)競爭等).其收斂速度與樣本數(shù)有關(guān)，與待求問題維度無關(guān)，而普通積分方法計算復(fù)雜度隨著維數(shù)的增加指數(shù)增加.因此對于很多高維計算，蒙特卡羅模擬是一個非常好的計算方式.

為達(dá)到較高的模擬精度和實(shí)現(xiàn)實(shí)時模擬，需要大量復(fù)雜計算并且計算次數(shù)可能達(dá)到上百萬次，因此需要對蒙特卡羅模擬進(jìn)行硬件加速.

現(xiàn)場可編程門陣列(FPGA)加速相對于通用CPU 具有高度并行、速度快、功耗低等優(yōu)點(diǎn);而相對于ASICs，則具有設(shè)計時間短、現(xiàn)場可編程等優(yōu)點(diǎn).因此FPGA 可重構(gòu)計算在蒙特卡羅模擬的硬件加速中得到了越來越多的應(yīng)用.

隨機(jī)數(shù)是蒙特卡羅模擬的重要模塊，可分為均勻分布隨機(jī)數(shù)和高斯分布隨機(jī)數(shù).均勻分布隨機(jī)數(shù)常用于高斯分布隨機(jī)數(shù)的產(chǎn)生，提供“隨機(jī)性”;高斯分布隨機(jī)數(shù)是應(yīng)用最多的分布，可用來近似很多隨機(jī)過程.高斯分布隨機(jī)數(shù)產(chǎn)生的方法有:反函數(shù)法、變換法、拒絕-接受法和遞歸法等［6-9］.反函數(shù)法和變換法精度較高，但計算復(fù)雜度高，硬件消耗大，設(shè)計復(fù)雜.拒絕-接受法和遞歸法具有結(jié)構(gòu)簡單、速率快的特點(diǎn)，但遞歸法輸出存在相關(guān)性，拒絕-接受法輸出速率不為常數(shù).

Marsaglia 等［10］提出的基于拒絕- 接受法的Ziggurat 方法具有結(jié)構(gòu)簡單、設(shè)計周期短的優(yōu)點(diǎn).文獻(xiàn)［11］對Ziggurat 方法進(jìn)行改進(jìn)后只保留了矩形區(qū)域，電路結(jié)構(gòu)簡單，克服了輸出不恒定的缺點(diǎn)，但對概率密度函數(shù)(PDF)分割精度不夠.文中在文獻(xiàn)［11］的基礎(chǔ)上，針對PDF 的尾部和頂部區(qū)域，采用嵌套分割的方法提高對尾部和頂部區(qū)域的近似精度，相對于傳統(tǒng)Ziggurat 方法減少了設(shè)計周期和硬件資源消耗.

1 高斯分布隨機(jī)數(shù)算法的改進(jìn)

1.1 Ziggurat 方法基本原理

Ziggurat 方法屬于拒絕-接受法，其基本原理如圖1 所示［12］.曲線p(x)=exp(-x2/2)，x ＞0 為高斯分布的概率密度函數(shù)，C 為曲線p(x)與坐標(biāo)軸構(gòu)成的區(qū)域.Z 由圖中的n-1 個矩形和1 個底部區(qū)域構(gòu)成，底部區(qū)域包括底部條形區(qū)域和尾部區(qū)域，尾部區(qū)域向右延伸至無限遠(yuǎn).其中，x(2)至x(n)為n-1個矩形的右側(cè)邊界系數(shù)，n 為構(gòu)成Z 的區(qū)域的個數(shù)，各子區(qū)域的面積均為v.

圖1 Ziggurat 算法示意圖Fig.1 Schematic diagram of Ziggurat algorithm

令i 為滿足條件1≤i≤n 的整數(shù)，當(dāng)i ＞1 時，為矩形區(qū)域Ri，當(dāng)i=1 時為底部區(qū)域.其原理為:

(1)當(dāng)i ＞1 時，各矩形區(qū)域Ri的右側(cè)邊界的坐標(biāo)為x(i)，Ri中一個隨機(jī)點(diǎn)的橫坐標(biāo)為x =U1x(i).若x ＜x(i-1)，即圖1 中子矩形區(qū)域，則其必然在C中，返回對應(yīng)的橫坐標(biāo)x.

(2)若x≥x(i-1)，且在曲線p(x)下方，滿足［p(x(i-1))-p(x(i))］U2＜p(x)-p(x(i))，則(x，y)來自契形區(qū)域，返回x 值.

(3)若i=1，且U3＜rf(r)/v，則認(rèn)為x 來自底部的條形區(qū)域，否則認(rèn)為來自尾部區(qū)域.

其中，U1、U2、U3為(0，1)上相互獨(dú)立的均勻分布隨機(jī)數(shù).

令區(qū)域Z 最右側(cè)矩形區(qū)域的邊界系數(shù)x(n)為r(即矩形區(qū)域最大輸出值)，文獻(xiàn)［11］計算得到:n取1 024 時，矩形區(qū)域占99.84%，契形區(qū)域占0.082%，尾部區(qū)域占0.086%.契形區(qū)域和尾部區(qū)域?qū)?yīng)的概率很小，卻需要計算指數(shù)、乘法、除法等復(fù)雜運(yùn)算，因此文獻(xiàn)［11］對Ziggurat 算法進(jìn)行簡化，只保留矩形區(qū)域，并將原算法的256 個分塊提高到1024 個，簡化了硬件結(jié)構(gòu).其不足是對契形區(qū)域的近似仍然不夠，并且輸出最大值r 僅為4.29.而在蒙特卡羅模擬中，雖然偏差較大的值出現(xiàn)概率很小，但對最終模擬結(jié)果影響較大.

1.2 改進(jìn)后的嵌套分割算法

文中對文獻(xiàn)［11］的方法進(jìn)行了改進(jìn):

(1)首先將待求函數(shù)的概率密度函數(shù)分割為若干個矩形區(qū)域.文中先將高斯分布概率密度函數(shù)分割為頂部、中部和尾部區(qū)域，再用嵌套分割的方法將頂部和尾部區(qū)域分割為矩形區(qū)域，圖2 為其分割示意圖，圖中x1(i)為x1層矩形邊界系數(shù).

圖2 嵌套分割示意圖Fig.2 Schematic diagram of nested segmentation

(2)產(chǎn)生概率密度函數(shù)與對應(yīng)子區(qū)域相同的隨機(jī)數(shù).文中采用Ziggurat 算法中矩形系數(shù)x(i)與(0，1)上均勻分布隨機(jī)數(shù)U 相乘的方法產(chǎn)生概率密度函數(shù)為矩形的隨機(jī)數(shù).

(3)選擇對應(yīng)子區(qū)域的隨機(jī)數(shù)作為輸出，其被選擇的概率正比于子區(qū)域的面積.

新算法的特點(diǎn)為:①省去了Ziggurat 算法中的拒絕-接受過程，輸出為常速率.②對于中部區(qū)域，將分區(qū)數(shù)從1024 提高到4096.實(shí)驗發(fā)現(xiàn)，分塊數(shù)的進(jìn)一步提高對頂部和尾部精度提高很小，因此需要單獨(dú)優(yōu)化.③對于尾部區(qū)域，采用嵌套分割方法，最大輸出值從4.29 提高到8.72.④對于頂部區(qū)域，由于每個塊的面積相等，當(dāng)系數(shù)絕對值較小時，矩形較高，存在較大誤差.采用嵌套分割算法改進(jìn)后系數(shù)差不到原有系數(shù)差的2%.

1.2.1 對頂部區(qū)域的近似

將文獻(xiàn)［11］中的分割數(shù)提高到4 096 后，其矩形區(qū)域邊界系數(shù)差lg［x(i)-x(i-1)］如圖3 所示，橫坐標(biāo)為矩形區(qū)域邊界系數(shù)序號i，縱坐標(biāo)為相鄰兩系數(shù)的對數(shù)差lg［x(i)-x(i-1)］.定義第一層為x1層，系數(shù)序號為［1，4096］，示意圖如圖2 所示，橫坐標(biāo)x1(1)至x1(4096)為x1層矩形區(qū)域的右側(cè)邊界系數(shù)，縱坐標(biāo)為概率密度函數(shù).由于矩形面積恒定，系數(shù)越小，矩形區(qū)域高度越大，誤差也越大.當(dāng)i ＜220 時，系數(shù)誤差大于10-4.因此對前256 個子區(qū)域進(jìn)行嵌套分割.將x1層前256 個子區(qū)域分割為1024個n1層子區(qū)域，每個子區(qū)域的大小為v×256/1024.

圖3 頂部區(qū)域前300 個矩形系數(shù)差Fig.3 Coefficient difference of the first 300 rectangles of top area

嵌套的第1 級為n1層，系數(shù)序號為［1，1 024］，其嵌套分割示意圖如圖4 所示，滿足x1(257)=n1(1025)，n1(1)=0.其中x1(257)為圖2 中x1層第257 個矩形右側(cè)邊界系數(shù)，n1(1 025)為圖4 中n1層第1025 個邊界系數(shù).對n1層前128 個子區(qū)域再次進(jìn)行嵌套分割得到1024 個n2層子區(qū)域.以此類推，不斷對系數(shù)進(jìn)行嵌套分割，使系數(shù)之差接近10-4.最終所得頂部區(qū)域只讀內(nèi)存ROM 地址如圖5 所示.

圖4 頂部區(qū)域嵌套分割示意圖Fig.4 Schematic diagram of nested segmentation of top region

圖5 頂部區(qū)域ROM 地址Fig.5 ROM address of top region

可通過以下公式用二分法計算邊界系數(shù)x(i):

尾部面積為

矩形面積為

解得

式中，r=x1(4096).

在迭代第一步，將區(qū)間［a，b］中點(diǎn)(a +b)/2 賦給r(a、b 為預(yù)估的r 的區(qū)間)，通過式(1)計算得到v，通過式(3)迭代計算x(i).當(dāng)i 在x1層迭代到i =257時，跳轉(zhuǎn)到n1層，并將v 替代為vn1=v ×256/1 024，其中vn1為n1層子區(qū)域面積.跳轉(zhuǎn)到n1層后再次將i 從1024 迭代到129，跳轉(zhuǎn)到n2層.以此類推，得到n4層第一個系數(shù)n4(1).對區(qū)間［a，b］采用二分法不斷迭代直到n4(1)=0.此時的r 即為區(qū)域Z 最右側(cè)第4096 個矩形邊界系數(shù)，即r=x1(4096)=4.39.

1.2.2 對尾部區(qū)域的近似

為提高輸出的高斯分布隨機(jī)數(shù)的最大值，需要對尾部區(qū)域進(jìn)行嵌套分割.如圖6 所示，將尾部區(qū)域分割為1 024 個子區(qū)域，通過式(3)從1 024 迭代到2，使得x2層第2 個系數(shù)等于x1層第4 096 個系數(shù)，即x2(2)=x1(4096)，得到x2層最大輸出值r2=x2(1024)=5.80.為進(jìn)一步提高最大輸出值r，繼續(xù)對x2層尾部進(jìn)行嵌套分割得到r3= x3(1024)=6.92，r4=x4(1024)=7.87，r5=x5(1024)=8.72….當(dāng)?shù)降? 層時，r 已經(jīng)超過8，因此文中采用四級迭代r=r5=8.72.相比于原算法的r=4.29，文中改進(jìn)后的r 提高了107.9%.最終所得尾部區(qū)域ROM 地址分層結(jié)構(gòu)如圖7 所示.

圖6 尾部區(qū)域嵌套分割示意圖Fig.6 Schematic diagram of nested segmentation of tail region

圖7 尾部區(qū)域ROM 地址Fig.7 ROM address of tail region

2 算法的硬件實(shí)現(xiàn)

新算法硬件系統(tǒng)結(jié)構(gòu)見圖8.

圖8 改進(jìn)算法系統(tǒng)框圖Fig.8 Block diagram of the proposed method

乘法器共有兩路輸入，α 輸入為16 位的(-1，1)上均勻分布隨機(jī)數(shù)，β 輸入為取自ROM 的16 位系數(shù)，兩者相乘得到32 位高斯分布隨機(jī)數(shù)，小數(shù)位為27 位.文中采用文獻(xiàn)［13］中Combined Tausworthe 方法產(chǎn)生均勻分布隨機(jī)數(shù).

2.1 多選器選擇信號的生成

圖8 中的ROM 選擇模塊(rom_select)如圖9所示.

圖9 ROM 選擇模塊結(jié)構(gòu)Fig.9 Structure of rom_select module

其中Tail_area 為尾部分割模塊，輸出sel_x2-sel_x5 分別是x2-x5層控制信號，rom_coef_1-rom_coef_5 分別為x1-x5層輸出的系數(shù)值;Top_area 為頂部區(qū)域模塊，各信號含義與Tail_area類似.Top_area 和Tail_area 輸出控制信號通過Concat 模塊組合為一個9 位信號，再通過ROM2 轉(zhuǎn)化為多路選擇器coef_mux 控制信號，選擇對應(yīng)的系數(shù).

多選器選擇信號通過ROM2 查找表生成，根據(jù)改進(jìn)后的算法，其對應(yīng)的真值見表1，0 代表低電平，1 代表高電平，X 代表0 或者1.

表1 ROM2 查找表Table 1 Look-up table of ROM2

2.2 ROM 選擇算法的實(shí)現(xiàn)

由改進(jìn)的算法可知，為選擇對應(yīng)ROM，需要對ROM 地址數(shù)據(jù)進(jìn)行判斷，以決定是尾部、頂部或者中部區(qū)域.下面以頂部區(qū)域為例介紹對應(yīng)算法的硬件實(shí)現(xiàn).

如圖10 所示為頂部區(qū)域模塊，自頂向下依次為n1-n4層，每層包括選擇信號生成模塊和系數(shù)存儲單元兩部分.以n1、n2層為例，其原理可描述如下.

圖10 Top_area 模塊示意圖Fig.10 Schematic diagram of Top_area module

(1)sel_logic_n1:選擇信號生成模塊，Slice2 提取x1層高4 位通過sel_logic_n1 進(jìn)行或非，生成n1層選擇信號sel_n1，sel_n1 高電平時(此時12 位x1層地址小于256，即高4 位為低電平)跳轉(zhuǎn)到n1層.

(2)n1層ROM:系數(shù)存儲單元，Slice1 提取輸入的n1層10 位地址信號輸出到“n1層ROM”中，取出對應(yīng)的系數(shù)，輸出到rom_coef_n1.

(3)n2層的Slice3 提取Slice1 中的高3 位進(jìn)行異或，生成n2層選擇信號sel_n2，sel_n2 低電平時輸出n1層系數(shù)rom_coef_n1，高電平時跳轉(zhuǎn)到n2層.

2.3 整體仿真結(jié)果

在System Generator 中對硬件實(shí)現(xiàn)后的整體電路進(jìn)行仿真，仿真結(jié)果如圖11 所示.其中urn 為乘法器輸入的均勻分布隨機(jī)數(shù)，coef 為rom 選擇模塊輸出的系數(shù)值，out 為乘法器輸出，即文中算法生成的高斯分布隨機(jī)數(shù).

圖11 電路仿真結(jié)果Fig.11 Simulation results of the circuit

sel_x2 為x2層控制信號，高電平時選擇尾部區(qū)域;sel_n1 為n1層選擇信號，高電平時選擇頂部區(qū)域;sel_n2 為n2層選擇信號.由表1 可知，3 個信號的組合“000”、“001”、“010”分別輸出“0”、“0”、“5”;當(dāng)rom_out 為“0”時，coef_mux 輸出x1層系數(shù)coef_x1;當(dāng)為“5”時，輸出頂部n1層系數(shù)coef_n1.表2為改進(jìn)后算法在Xilinx Virtex 4 xc4vls100 FPGA 上的綜合結(jié)果，速率為270 MHz，最大輸出值為8.72.

文中的嵌套分割隨機(jī)數(shù)發(fā)生器所消耗的硬件資源比文獻(xiàn)［8］基于Ziggurat 算法的硬件結(jié)構(gòu)減少四分之三.在并行蒙特卡羅模擬中，若將一個蒙特卡羅模擬在M 個蒙特卡羅模擬加速器中加速，則其計算時間為單個加速器的1/M.每個加速器中隨機(jī)數(shù)均為完全相同的分布，因此RAM 中的系數(shù)值也完全相同.而FPGA 中雙口RAM 可同時讀出兩路數(shù)據(jù)，因此一組RAM 可以供兩路隨機(jī)數(shù)發(fā)生器使用，平均每路隨機(jī)數(shù)發(fā)生器消耗12/2 =6 個塊RAM，以“6(12)”列于表2 最后一行.

本算法的特點(diǎn)為:不僅可以生成高斯分布隨機(jī)數(shù)，還可以生成概率函數(shù)從最高點(diǎn)向兩邊遞減的分布函數(shù);輸出隨機(jī)數(shù)最大值可以通過增加尾部嵌套層數(shù)實(shí)現(xiàn)，文中通過四級嵌套最大輸出可以達(dá)到8.72.

表2 綜合結(jié)果對比Table 2 Comparison of synthesis results

3 輸出結(jié)果的統(tǒng)計檢驗

3.1 Diehard 檢驗

Diehard 檢驗包含有一系列檢驗［17］，包括游程檢驗、DNA 檢驗、生日檢驗等.檢驗時選擇檢驗的種類，程序返回一個P 值，如果輸入文件包含的隨機(jī)數(shù)為獨(dú)立的隨機(jī)數(shù)，則P 值應(yīng)該滿足在［0，1)上均勻分布.

假設(shè)x、y 服從正態(tài)分布，則

為(0，1)上均勻分布隨機(jī)數(shù).將文中嵌套分割算法所生成的600 萬個高斯分布隨機(jī)數(shù)轉(zhuǎn)換為300 萬個均勻分布隨機(jī)數(shù)，部分檢驗結(jié)果見表3，所生成的均勻隨機(jī)數(shù)通過了Diehard 檢驗.

表3 Diehard 檢驗結(jié)果Table 3 Results of Diehard tests

3.2 χ2擬合檢驗

χ2檢驗可以檢驗樣本的概率密度函數(shù)與理論函數(shù)的差異，檢驗統(tǒng)計量為

式中，Oi和Ei分別為區(qū)間i 上實(shí)際觀察和理論計算得到的樣本數(shù)，Q 為區(qū)間個數(shù).原假設(shè)為樣本x 服從正態(tài)分布，在原假設(shè)成立的條件下，計算ξ2對應(yīng)的P值，若P 小于置信度(一般為0.05)，則拒絕原假設(shè)，否則接受原假設(shè).每次測試100 萬個隨機(jī)數(shù)，5 次結(jié)果見表4，均通過了統(tǒng)計檢驗.

表4 卡方檢驗結(jié)果Table 4 Results of Chi-squared test

其中“Matlab”為Matlab 自帶“randn”函數(shù)所產(chǎn)生的隨機(jī)數(shù)

3.3 K-S 檢驗

K-S 測試可以檢驗樣本的實(shí)際積累分布函數(shù)與理論分布函數(shù)的差異［18］，檢驗統(tǒng)計量為max(|F(x)-G(x)|).

原假設(shè)為樣本x 服從正態(tài)分布，若原假設(shè)為真，則其值應(yīng)該較小，否則應(yīng)該拒絕.每次測試5 萬個隨機(jī)數(shù)，5 次結(jié)果見表5，通過了K-S 檢驗.

表5 K-S 測試結(jié)果Table 5 Results of K-S test

3.4 相關(guān)性檢驗

如果相鄰輸出存在相關(guān)性，則其輸出存在規(guī)則的晶格結(jié)構(gòu)［16］.文中算法產(chǎn)生的隨機(jī)數(shù)散點(diǎn)圖(Ui，Ui+1)如圖12 所示，沒有明顯的相關(guān)性.

在Matlab 中計算改進(jìn)后算法生成的兩萬個隨機(jī)數(shù)的±2000 個lag 對應(yīng)的自相關(guān)函數(shù).如圖13 所示，除原點(diǎn)外其他點(diǎn)歸一化后的自相關(guān)函數(shù)接近0，沒有明顯的相關(guān)性.

圖12 生成隨機(jī)數(shù)的散點(diǎn)圖Fig.12 Scatter plot of generated random numbers

圖13 生成隨機(jī)數(shù)的自相關(guān)函數(shù)Fig.13 Autocorrelation function of generated random numbers

4 結(jié)語

文中針對Ziggurat 算法輸出速率不恒定、算法中出現(xiàn)概率很小的支路卻占用大量的硬件資源和設(shè)計復(fù)雜等問題，對高斯分布的概率密度函數(shù)進(jìn)行矩形嵌套分層分割，對每一個矩形塊采用Ziggurat 中算法生成概率密度函數(shù)為矩形的隨機(jī)數(shù).針對仿真中出現(xiàn)的極值情況，對尾部區(qū)域進(jìn)行了優(yōu)化處理，最大輸出可以達(dá)到8.72.統(tǒng)計檢驗結(jié)果表明，新嵌套分割算法生成的隨機(jī)數(shù)為獨(dú)立不相關(guān)的高斯分布隨機(jī)數(shù).

［1］Luu J，Redmond K，Lo W C Y，et al.FPGA-based Monte Carlo computation of light absorption for photodynamic cancer therapy［C］∥Proceedings of 17th IEEE Symposium on Field Programmable Custom Computing Machines.Napa:IEEE Computer Society，2009:157-164.

［2］Simon A，Tudoran R，Cret O，et al.FPGA-based Monte-Carlo computation of the electric potential in homogeneous and non-homogeneous spaces［C］∥Proceedings of the 8th FPGA World Conference.New York:ACM，2011:63-68.

［3］Castillo J V，Atoche A C，Longoria-Gandara O，et al.An efficient Gaussian random number architecture for MIMO channel emulators［C］∥Proceedings of IEEE Workshop on Signal Processing Systems.Beirut:IEEE，2011:316-321.

［4］Sanchez-Roman D，Moreno V，Lopez-Buedo S，et al.FPGA acceleration using high-level languages of a Monte-Carlo method for pricing complex options［J］.Journal of Systems Architecture，2013，59(3):135-143.

［5］de Schryver C，Shcherbakov I，Kienle F，et al.An energy efficient FPGA accelerator for Monte Carlo option pricing with the Heston model［C］∥Proceedings of International Conference on Reconfigurable Computing and FPGAs.Cancun:IEEE，2011:468-474.

［6］Lee D U，Cheung R C C，Villasenor J D，et al.Inversionbased hardware Gaussian random number generator:a case study of function evaluation via hierarchical segmentation［C］∥Proceedings of IEEE International Conference on Field Programmable Technology.Bangkok:IEEE，2006:33-40.

［7］Fung E，Leung K，Parimi N，et al.ASIC implementation of a high speed WGNG for communication channel emulation［C］∥Proceedings of IEEE Workshop on Signal Processing Systems Design and Implementation.Austin:IEEE，2004:304-309.

［8］Zhang G，Leong P H W，Lee D U，et al.Ziggurat-based hardware Gaussian random number generator［C］∥Proceedings of International Conference on Field Programmable Logic and Applications.Tampere:IEEE，2005:275-280.

［9］Lee D U，Luk W，Villasenor J D，et al.A hardware Gaussian noise generator using the Wallace method［J］.IEEE Transactions on Very Large Scale Integration (VLSI)Systems，2005，13(8):911-920.

［10］Marsaglia G，Tsang W W.The Ziggurat method for generating random variables［J］.Journal of Statistical Software，2000，5(8):1-7.

［11］王靈芝.用FPGA 產(chǎn)生高斯隨機(jī)數(shù)及其在量子密碼中的應(yīng)用［D］.福州:福州大學(xué)物理與信息學(xué)院，2006:17-19.

［12］Thomas D B，Luk W，Leong P H W，et al.Gaussian random number generators ［J］.ACM Computing Surveys(CSUR)，2007，39(4):11/1-38.

［13］L’Ecuyer P.Maximally equidistributed combined Tausworthe generators ［J］.Mathematics of Computation，1996，65(213):203-213.

［14］Lee D U，Villasenor J D，Luk W，et al.A hardware Gaussian noise generator using the Box-Muller method and its error analysis［J］.IEEE Transactions on Computers，2006，55(6):659-671.

［15］Thomas D B，Luk W.Efficient hardware generation of random variates with arbitrary distributions［C］∥Proceedings of IEEE Symposium on Field-Programmable Custom Computing Machines.Napa:IEEE，2006:57-66.

［16］Malik J S，Malik J N，Hemani A，et al.Generating high tail accuracy Gaussian random numbers in hardware using central limit theorem［C］∥Proceedings of IEEE/IFIP 19th International Conference on VLSI and Systemon-Chip.Hong Kong:IEEE，2011:60-65.

［17］L' Ecuyer P.Random number generation［M］.Heidelberg:Springer，2012:21-23.

［18］Knuth D E.The art of computer programming，vol.2:seminumerical algorithms［M］.Massachusetts:Addison_Wesley，1981:48-55.