王 超,張秋艷,張 姍,王 龍
(中國電子信息產業集團有限公司第六研究所,北京 100083)
隨機數廣泛應用于概率統計、模擬仿真、信息安全和數字通信等諸多領域。當需要產生海量的隨機數時,傳統串行的隨機數生成算法(隨機數發生器,Random Number Generator)時間過長,難以達到實際需求。
本文在基于線性反饋移位寄存器(Linear Feedback Shift Register,LFSR)產生偽隨機數的基礎上,利用采樣定理,提出了一種基于多核處理器的新算法。在新算法中,將串行產生方式改為并行產生方式,提高了產生偽隨機數的效率,并且新算法具有并行與串行結果一致的特性,即新算法保持了通用性。本文首先證明了新算法的可計算性、確定性和結果一致性,然后給出了軟件實現流程和硬件推廣分析,最后在Intel(R) Core(TM) 四核CPU i7-6700上進行偽隨機數生成實驗,相對于傳統的串行算法,加速比已經接近4。
并行計算機[1]是相對于串行計算機而言的,所謂串行計算機(順序計算機)就是單個處理單元順序執行計算機程序的計算機。典型的并行計算機有多核處理器、現場可編程門陣列(Field-Programmable Gate Array,FPGA)芯片和圖形處理器(Graphics Processor Unit,GPU)。
多核處理器,又稱為單芯片多處理器(Chip Multiprocessors,CMP),其各個處理器并行執行不同的任務,通過線程并行性來取代越來越復雜的指令集并行,以此提高處理器的性能。
FPGA芯片是一種擁有高密度數字電路、高處理性能和可編程使用的信號處理器件,其通過消耗內部邏輯資源塊實現并行處理。
以CUDA(Computer Unified Device Architecture)平臺為代表的圖形處理器,具有相當高的內存帶寬以及大量的浮點計算單元,其通過使用大量線程來充分利用多計算核心,從而實現高性能。
對于并行計算機,盡管它能夠提高多任務系統的性能,但是它不能提升串行系統的性能。因此,如果現有串行算法設計思想不加以改變,將無法充分利用并行計算機的計算能力。
隨機數包括物理隨機數(真隨機數)和偽隨機數兩類,若不特別說明,本文所涉及的隨機數都是指偽隨機數。
隨機數[2]可按均勻性劃分為均勻隨機數和非均勻隨機數。均勻隨機數是產生非均勻隨機數的基礎,如正態分布、指數分布等隨機數都可以用均勻隨機數經過變換得到。當今流行的均勻隨機數序列有:反饋移位寄存器(m序列、M序列)、二次剩余序列、霍爾(Hall)序列、孿生素數序列、混沌映射序列、進位加法和借位減法序列[3-4]等。其中反饋移位寄存器包括線性遞歸序列和非線性遞歸序列兩類,由于非線性遞歸序列的周期特性和統計特性還沒有成熟結論,分析這類序列密碼的安全性比較困難,因而LFSR是序列密碼設計中常用的一種初始亂源。
通過將種子密鑰作為LFSR的初態,按照遞推關系,產生一個周期長、統計特性好的初始亂源,從而為進一步利用非線性函數、有記憶變換等設計手段,產生抗破譯能力強的隨機數序列提供“原料”。
2.1.1定義與定理
為證明本文的新隨機數生成算法(簡稱新算法)具有可計算性和確定性,以及并行與串行結果的一致性[5],下面先引入幾個定義和定理,然后通過具體的算法過程來證明。

定義2:如果二元域上的n級線性遞歸序列的周期是2n-1,則稱該序列是二元域上的一條n級m序列,并稱其對應的反饋多項式是本原多項式。


定理1:二元域上的n級m序列的游程特性:在一個周期內,不存在長度大于n-1的0游程。
定理2:n級m序列的2i采樣(i=0,1,…,n-1)都與序列平移等價。
定理3:n級線性反饋移位寄存器由n個連續項(初始狀態)和反饋多項式完全確定。
定理4:n次本原多項式的狀態圖由1個長度為1的圈(由零狀態構成)和1個長度為2n-1的圈組成。
2.1.2算法具體過程
設f(x)是二元域上的n次本原多項式,S0是m序列的初始狀態(非全0向量,通常取值為(0,…,0,1)),并行線程個數記為2r。
步驟1:由f(x)和S0,通過基于乘法電路設計的線性反饋移位寄存器,生成長度為n×2r的輸出序列,記為序列a。

步驟3:將2r個線性反饋移位寄存器輸出進行拼接得到最終輸出。
2.1.3并行原理
(1)證明新算法具有可計算性
①步驟2中2r個初始狀態與步驟1中序列a的前n×2r項相同,如公式(1)所示:
(1)


綜上可知步驟2具有可計算性,于是新算法具有可計算性。
(2)證明新算法具有確定性
定理3保證2r個線性反饋移位寄存器都完全唯一確定,因此步驟2具有確定性,于是證明新算法具有確定性。
(3)證明新算法具有并行與串行結果一致性
序列a可以通過序列a,La,L2a,…,L2r-1a的2r采樣拼接而成,如公式(2)所示。
(2)
證畢。
2.2.1組成與功能
新算法由主控線程和工作線程兩部分組成,總體架構如圖1所示。

圖1 新算法總體架構圖
(1)初始化與預計算模塊
該模塊主要完成:①工作線程個數2r的設置;②以下變量的初始化:用于索引的全局鏈表及保護它的線程互斥量,待分配任務計數、已運行任務計數和啟動線程結束標志及保護三個變量的工作線程互斥量、工作線程條件變量;③反饋多項式的設置,預生成長度為n×2r的輸出序列,作為工作線程的輸入參數。
(2)工作線程創建模塊
該模塊完成創建2r個工作線程,其參數為反饋移位寄存器的初始狀態和工作線程序號。
(3)任務啟動與同步模塊
該模塊觸發工作線程遷移出等待階段,當工作線程完成任務時,判定是否當前批次任務都完成,若未完成則繼續等待。
(4)啟動工作線程結束模塊
該模塊觸發工作線程遷移出等待階段,允許工作線程結束。
(5)等待工作線程結束模塊
該模塊等待所有工作線程結束。
(6)銷毀與回收模塊
該模塊銷毀互斥量和條件變量,回收運行過程中申請的堆內存。
(7)工作線程模塊
該模塊由等待階段、領取任務同步、執行任務階段、完成任務同步和結束本線程組成。
2.2.2線程同步信息
圖1中以黑色實心箭頭形式標出主線程與工作線程之間的同步信息:主控線程中“任務啟動與同步模塊”觸發工作線程中“等待階段”;主控線程中“啟動工作線程結束模塊”觸發工作線程中“等待階段”;工作線程中“完成任務同步”觸發主控線程中“任務啟動與同步模塊”。圖1中以空心箭頭形式標出工作線程之間的同步信息:工作線程中領取任務同步觸發其他工作線程的等待階段。
2.2.3線程同步技術
考慮到各工作線程處理任務的相似性,得出各工作線程處理任務消耗的時間接近,于是減化同步控制,僅使用線程互斥量THREAD_MUTEX和工作線程條件THREAD_COND控制所有工作線程同步。本算法在實現中使用了最低數量的線程共享變量(待分配任務計數、已運行任務計數和啟動線程結束標志)。全局鏈表與線程共享變量不存在資源競爭,為提高各工作線程效率,沒有使用THREAD_MUTEX保護全局鏈表,而是增加互斥量gPPHEAD_MUTEX來保護全局鏈表。
任務啟動與同步模塊通過設置待“分配任務計數”,觸發各工作線程開始執行任務,然后直到“已運行任務計數”達到要求,再繼續執行。
當完成任務啟動與同步后,此時工作線程可能仍處于執行任務階段,主控線程通過設置“啟動線程結束標志”告知各工作線程后續沒有新任務,可以結束本工作線程。
各工作線程通過判斷“分配任務計數”決定是等待還是執行任務;領取任務后“待分配任務計數”減一,并通知其他工作線程的等待階斷;完成任務后“已運行任務計數”加一,并通知主控線程的任務啟動與同步模塊,主控線程判斷“已運行任務計數”達到要求,決定是否繼續等待。
為了實現隨機數生成算法的并行與串行結果一致,傳統串行算法的流程如下:
(1)隨機數生成算法由n個模塊組成;
(2)第i(i=1,…,n)個模塊模擬運行i輪LFSR;
(3)內部狀態步進n輪,轉至上一過程。
本文提出的新算法,對應的流程如下:
(1)隨機數生成算法由n=2r個模塊組成;
(2)第i(i=1,…,n)個模塊模擬運行1輪LFSR;
(3)內部狀態步進1輪,轉至上一過程。
在FPGA[6]上的傳統串行算法,當n較大時,其第n個模塊相較第1個模塊復雜很多,使得FPGA的時鐘頻率無法很高。而新算法中,所有n個模塊的復雜度一樣,使得FPGA的時鐘頻率可以相對更高。
接下來考查GPU[7-9],例如GeForce8800GTX擁有128個執行單元,為充分利用,需要多線程。同時,由于CUDA平臺中對全局內存的存取有較大延遲,線程會因為無法及時讀取或寫入數據而處于等待狀態,需通過使用超量的線程來隱藏全局內存延遲,從而需要一個較高的計算和存取比率[10]。于是,為充分利用GPU,需要數千個線程。傳統串行算法中所有n個模塊的算法邏輯不同,而新算法中所有n個模塊的算法邏輯是一樣的,僅初始值不同,并且初始值為少量輸入,GPU恰好適合新算法的加速實現。
綜上所述,本文提出的新算法思想不僅適用于多核處理器,也同樣適用于FPGA和GPU。
本文采用的測試平臺是Intel(R) Core(TM) 四核CPU i7-6700,主頻3.40 GHz,內存8 GB,搭載Windows 7(Service Pack 1)旗艦版操作系統和GCC 4.6編譯環境,編譯優化選項為O3。
本文采用
為了降低對時間的測量誤差,本文采用2 000作為單個測試的倍數。通過計算平均值獲取平均情況,采用8次單個測試作為1個測試組。為了消除內存使用量對算法性能的影響,運行1組測試后,立即進行堆內存釋放,然后再進行下一組測試,累計進行4組測試。對于1套給定參數,總計產生32個運行時間和1個平均值。
線性反饋移位寄存器的參數包括反饋函數和初始狀態,其中反饋函數的次數、反饋函數的非0項數和初始狀態都影響隨機數生成算法的速度。本文采用常見的168次反饋函數。初始狀態需為非全0向量,本文采用常見的(0,…,0,1)。為考查不同反饋函數項數對測試結果的影響,本文使用2套參數。第1套參數,取項數最少的x168+x17+x15+x2+1作為反饋函數。通過采樣定理(取采樣間隔為5,5為最小的互素數)和Barlekamp Massey算法得到另一個項數為39的本原多項式x168+x111+x110+x109+x98+x93+x82+x80+x71+x70+x69+x64+x63+x62+x58+x53+x51+x46+x43+x41+x40+x39+x36+x35+x34+x32+x30+x29+x28+x25+x24+x23+x22+x21+x18+x16+x15+x2+1作為第2套參數。
經過4組×8次/組測試,新算法并行相對串行加速比如圖2所示,當工作線程個數為1時,相對于傳統的串行算法,本文提出的新算法的加速比接近1;當工作線程個數為2時,加速比接近2;當工作線程個數為4或64時,加速比接近4。

圖2 新算法并行相對串行加速比圖
綜上所述,本文提出的新算法的加速比最大值與測試平臺CPU的物理核心個數一致,即在Intel(R) Core(TM) CPU i7-6700上新算法的加速比最大為4,此結論與理論推導一致。同時,驗證新算法對反饋函數的非0項個數不敏感,于是新算法適用于各種線性反饋移位寄存器參數。
本文主要設計和實現了一種基于LFSR的適用于多核處理器的新偽隨機數生成算法。新算法在并行運行時與經典串行算法產生一致的隨機數,相對于傳統的串行算法,加速比可以達到多核處理器的物理核心數,同時保持了通用性。下一步研究工作是將此算法向GPU與FPGA推廣。