Python下探究隨機數的產生方法
資源下載
#本文PDF版下載
Python下探究隨機數的產生方法(或者單擊我博客園右上角的github小標,找到lab102的W7目錄下即可)
#本文代碼下載
幾種隨機數算法集合(和下文出現過的相同)
前言
我們對於隨機數肯定不會陌生,隨機數早已成為了我們經常要用到的一個方法,比如用於密碼加密,數據生成,蒙特卡洛算法等等都需要隨機數的參與.那麽我們的電腦是怎麽才能夠產生隨機數的呢?是電腦自己的物理存在還是依靠算法?它到底是如何工作的呢?所以我也對這些問題有著好奇心,所以找到了許多資料學習了一下,現在就整理下將自己的理解分享給大家.
隨機數的概念
這裏我們先來看一下數學上對隨機數這個詞的定義,定義如下:
在連續型隨機變量的分布中,最簡單而且最基本的分布是單位均勻分布.由該分布抽取的簡單子樣稱為隨機數序列,其中每一個體稱為隨機數.單位均勻分布即[0,1]上的均勻分布.
由隨機數序列的定義可知,ξ1,ξ2,…是相互獨立且具有相同單位均勻分布的隨機數序列.也就是說,獨立性、均勻性是隨機數必備的兩個特點.
隨機數的真偽
為什麽說隨機數會分成真和偽呢?那是因為,真正的隨機數,是物理中量子力學的概念,而我們如果想要真正做到真隨機,那麽就需要真正的獨立於電腦的特殊設備來實現生成真正的隨機數(比如英特爾曾經的芯片組是通過用熱噪聲放大後,影響電壓控制的振蕩器並用另外一個高頻振蕩器來接受數據得到隨機數)以及往大了來說自然界的一切都可以用於作為隨機數發生器的因子.而另一種通過算法實現的被稱為偽隨機的方法是指根據算法和指定的不確定因素(也就是種子seed),我們在電腦中常會用到電腦時間來做為seed的值,但是這樣的話,除非是毫秒級甚至以下的單位,否則在很短的時間內生成大量的隨機數,很容易發生時間重合,造城最終隨機數相同的情況.所以,在偽隨機數中,種子(seed)的指定就顯得相當的重要了.據說UNIX系統中,將系統時間,連入WIFI,甚至按下的鍵盤次數都量化為了seed,參考指標越多,偽隨機數就越接近真正的隨機生成.
產生隨機數的方法
產生偽隨機數有很多的方法,而在本文中,將分別對其中的幾種偽隨機數生成方法進行介紹和實踐.方法如下:
- 線性同余方法
- 平方取中方法
- 梅森旋轉隨機生成法
- BOX-MULLER法
隨機數產生方法好壞的判定
對於這個問題,我在網上搜到了相當之多的判定方法,這之中有許多的判定標準和判斷工具.比如像NIST測試方法,TestU01的測試工具等等,這裏我就貼上德聯信息安全工作室的一套方法(其他方法作為參考內容會放在文末的文獻參考大家也可以去看看了解一下相關的知識)
德國聯邦安全辦公室(下簡稱德聯安全辦)總共公布了四個等級,這裏就稍微翻譯一下大致的意思:
- 包含相同隨機數序列的概率很低.
- 根據指定的統計檢驗區別於“真隨機數”的數字序列.
- 它不可能被任何人能夠通過計算得到,或以其他的方法進行猜測,得到任何給定的一個子序列,以及任何工作中產生的一切序列中的值,以及叠代器的工作狀態.
-
對於任何目的而言,任何人不能從隨機數生成器的內部狀態計算或猜測得到序列中的任何一個之前的數字/序列以及任何之前的隨機數的狀態.
任意範圍隨機數生成原理
我們經常會對隨機數有一定的要求範圍,那麽隨機數產生的一般步驟是怎樣的呢?因為真隨機涉及到了物理的量子.故本文全文只討論偽隨機數的生成方法,我們在python 中的random庫中,經常會用到如randint之類的方法來生成一定範圍內的隨機數.這之中主要用到的方法步驟為->指定一個特定的seed種子,根據seed再通過特定的隨機數產生方法,就可以做到在[0,1]這個範圍中取到隨機分布的隨機數,然後一定的變換,我們就可以得到一定範圍內的隨機數了.
偽隨機數的算法
O平方取中法生成隨機數
平方取中法的思路是:
- 選擇一個m位數N作為SEED種子->做平方運算(記(N* N))
- 判斷N*N的長度,若結果不足2*m個位,在最前面補0.
- 在這個N*N的數選中間m個位的數作為隨機數.
公式表示為:
隨機數為:
我們來看看平方取中的代碼實現:
#平方取中 -xlxw from time import time def rander(seed,n): if n ==1: return 0 seed = int(seed) length = len(str(seed)) seed = int(seed**2/pow(10,(length/2))) % int(pow(10.0,length)) print(str(seed) +" " ,end=‘‘) rander(seed,n-1) def main(): seed = time() rander(seed,100) main()
生成的隨機數的截圖:
平方取中法的缺點:
- 周期很短
- 分布並不均勻
但是它因為計算迅速所以也常被用於隨機數據的生成.用於試驗等地方.和它相近的方法還有:常數取中法和乘法取中法等
O線性同余法計算隨機數
先來解釋一下名字的含義,我們把線性同余拆開來理解.線性的意思是滿足可加性和齊次性的映射;同余的意思就是如果給定一個數M,使兩個整數a,b可以滿足(a-b)能被M整除,則稱a,b對M同余.另一種理解方法就是a和b同時被M除後取余數相等則為同余.
根據這種方法的兩個特性,我們可以作如下兩個式子:
(a +/- b) mod M = ((a mod M) +/- (b mod M)) mod m
a*b mod M = ((a mod M) * (b mod M) mod M
我們來看一下這個式子的組成部分的含義:
- 我們可以從這個式子中看出,這是一個遞歸公式,所以X0為它的基,X0就是我們要給出的Seed種子.X0又叫做初始值.
- A的定義是系數
- C的定義是是增數
- M的定義是模數,也可以看作循環周期(最大)
- C,M,A是影響生成隨機數質量的重要因素.
我們來看一下下面用Python寫的模擬線性同余方法的代碼:(由於每次只生成一個隨機數,所以沒有用到遞歸,seed用時間來表示)
這裏M/A/C的值用了glibc (used by GCC) 編譯器的參數值
M = 2^32 A = 1103515245 C = 12345
#Random - LCG Ver by xlxw #import from time import time,clock #define m = 2**32 a = 1103515245 c = 12345 def LCG(seed): seed = (a * seed + c) % m return seed/float(m-1) def main(): br = input("請輸入隨機數產生的範圍(用,隔開):") mi = eval(br.split(‘,‘)[0]) ma = eval(br.split(‘,‘)[1]) seed = time() rd = LCG(seed) ourd = int((ma-mi)*rd) + mi print("隨機生成的數字式:{}".format(ourd)) main()
程序運行截圖:
通過一個seed生成多個隨機數的代碼:
#Random - LCG Ver by xlxw #import from time import time,clock #define m = 2**32 a = 1103515245 c = 12345 rdls = [] def LCG(seed,mi,ma,n): if n == 1: return 0 else: seed = (a * seed + c) % m rdls.append(int((ma-mi)*seed/float(m-1)) + mi) LCG(seed,mi,ma,n-1) def main(): br = input("請輸入隨機數產生的範圍(用,隔開):") co = eval(input("請輸入需要產生的隨機數的個數:")) mi = eval(br.split(‘,‘)[0]) ma = eval(br.split(‘,‘)[1]) seed = time() LCG(seed,mi,ma,co) print("隨機生成的數字",rdls) main()
程序運行截圖:
上面我們已經可以用線性同余來得到隨機數了.那麽我們來通過下面來看看生成的隨機數的概率是不是相同的,我們這裏順便測試一下生成隨機數的速度.
請看下面代碼:
#Random - LCG Ver by xlxw #import from time import time,clock import sys sys.setrecursionlimit(20000) #define m = 2**32 a = 1103515245 c = 12345 rdls = [] def LCG(seed,mi,ma,n): if n == 1: return 0 else: seed = (a * seed + c) % m rdls.append(int((ma-mi)*seed/float(m-1)) + mi) n = n-1 LCG(seed,mi,ma,n) def POSCheck():#得到各個數字出現次數 counts = {} for i in rdls: if i in counts: counts[i] = counts.get(i,0) + 1 else: counts[i] = 1 print(counts) def main(): br = input("請輸入隨機數產生的範圍(用,隔開):") co = eval(input("請輸入需要產生的隨機數的個數:")) mi = eval(br.split(‘,‘)[0]) ma = eval(br.split(‘,‘)[1]) seed = time() LCG(seed,mi,ma,co) POSCheck() #print("隨機生成的數字",rdls) main()
程序運行截圖:
(由於IDLE的遞歸深度限制在了1000以內,我通過修改了這個限制使得能遞歸20000次,再高則會失去響應,可能有什麽限制,所以這裏就拿生成20-29的隨機數20000運行三次來大致看看是否均勻)
我們可以看到出現次數大致上都分布在2000次左右,所以差不多可以看作產生的隨機數的概率是差不多一樣的.
M/A/C的選定(為了生成高質量隨機數)
- M:對於隨機數的生成而言,隨機數序列周期越長, 它就能夠具有在[0, 1] 上均勻分布及相互獨立.M越大周期就越大. 所以我們使得M接近於計算機能表示的最大的整數, 所以我們選擇2^32作為M的值
- A:對於乘數而言,1.對於M的任何一個質因子P,a-1能被P整除.2.如果4是M的因子,則a除以4余1(資料中提到,證明在數論中有)
-
C:增量C和M互質(證明同樣在數論中)
以上條件可以用於生成高質量隨機數,也就是將隨機數的循環達到M的值(其他情況下都是小於循環M)
線性同余法的缺點
- 線性同余法生成的隨機數周期太短.
- 如果被知道一定長度的序列就能破解出所有的隨機數生成序列(梅森旋轉也同樣有這種缺點)
O梅森旋轉法
梅森旋轉算法是一個偽隨機數發生算法.由松本真和西村拓士開發,是一個基於有限二進制字段上的矩陣線性遞歸的算法.可以快速產生高質量的偽隨機數,修正了樸素隨機數生成的缺陷.
梅森旋轉法是通過線性反饋移位寄存器來生成隨機數的.線性反饋移位寄存器-LFSR,是指給定前一狀態的輸出,將該輸出的線性函數再用作輸入的移位寄存器也就是對寄存器的某些位進行異或操作後作為輸入,再對寄存器中的各比特進行整體移位.而LFSR由兩個部分組成:
- 級數-一級一比特
- 反饋函數
解釋:
- 一:異或操作
異或運算通常用符號"⊕"表示,其運算規則為:
0⊕0=0 0同0異或,結果為0
0⊕1=1 0同1異或,結果為1
1⊕0=1 1同0異或,結果為1
1⊕1=0 1同1異或,結果為0
- 二:線性反饋移位寄存器的圖示:
- 三:某些位的說明
在LFSR線性反饋移位寄存器的解釋中說的“某些位”進行異或處理是由特征多項式來決定的.
- 四:特征多項式的要求
一個n階的本原多項式,在LFSR中要求使用本原多項式的原因是本原多項式才能使線性反饋移存器的周期達到最大.而本原多項式指的是:一個n次不可約多項式,如果只能整除1+x^2^n-1而 不能整除其它1+x^L(L<2^n-1),則這種不可約多項式就稱為本原多項式.
線性反饋移位寄存器的工作原理
例.以一個4位的線性反饋移位寄存器為例說明它的工作原理
反饋函數為:,可根據反饋函數的多項式作出線性移位寄存器邏輯框圖(如下圖):
設輸入的初始狀態為:
那麽計算步驟為:
(範例)
A3 A2 A1 A0
1 0 0 0
進行第一次異或判斷
進行移位處理,並將a3’的狀態給a3
A3 A2 A1 A0
1 1 0 0
之後再重復進行這個操作產生周期序列:
A3 A2 A1 A0
1 0 0 0
1 1 0 0
1 1 1 0
1 1 1 1
到這時,再進行異或判斷時
所以接下去的序列為:
A3 A2 A1 A0
1 0 0 0
1 1 0 0
1 1 1 0
1 1 1 1
0 1 1 1
再繼續進行:
1 0 1 1
0 1 0 1
1 0 1 0
1 1 0 1
0 1 1 0
0 0 1 1
1 0 0 1
0 1 0 0
0 0 1 0
0 0 0 1
1 0 0 0
我們可發現最後一行的序列和我們的初始狀態序列是一樣的.這就是一次移位寄存器的完整周期.可以看出這個序列的周期為15.在這一個周期裏有1-15所有整數,而且並沒有以固定的順序出現,說明了有很好的隨機性.
梅森素數Python的實現
#梅森旋轉算法 -xlxw #參考:mersenne twister from wikipedia #import from time import time import numpy as np #var index = 624 MT = [0]*index # MT[0] ->seed def inter(t): return(0xFFFFFFFF & t) #取最後32位->t def twister(): global index for i in range(624): y = inter((MT[i] & 0x80000000) +(MT[(i + 1) % 624] & 0x7fffffff)) MT[i] = MT[(i + 397) % 624] ^ y >> 1 if y % 2 != 0: MT[i] = MT[i] ^ 0x9908b0df index = 0 def exnum(): global index if index >= 624: twister() y = MT[index] y = y ^ y >> 11 y = y ^ y << 7 & 2636928640 y = y ^ y << 15 & 4022730752 y = y ^ y >> 18 index = index + 1 return inter(y) def mainset(seed): MT[0] = seed #seed for i in range(1,624): MT[i] = inter(1812433253 * (MT[i - 1] ^ MT[i - 1] >> 30) + i) return exnum() def main(): br = input("請輸入隨機數產生的範圍(用,隔開):") mi = eval(br.split(‘,‘)[0]) ma = eval(br.split(‘,‘)[1]) so = mainset(int(time())) / (2**32-1) rd = mi + int((ma-mi)*so) print("產生的隨機整數為:",rd) main()
程序運行截圖:
另一類偽隨機數生成方法
這一類隨機數還是偽隨機數.但是和上面幾種算法生成的偽隨機數不同在於上面幾種方法生成的偽隨機數追求的都是分布均勻.而下面講的是用於一些特殊情況下隨機數的生成方法:先讓我們來看看以下幾個例子:
- 當我們要對一個班的成績進行模擬實驗的時候,最貼近實際的是一種接近正態分布的成績分布情況.並且數據要隨機.
- 還有像模擬大學生體重情況,也和上面的一樣遵循類正態分布.
-
還有像產品的良品率,降雨量等等都是
所以在現實中我們需要符合正態分布的隨機數.而下面我們對其中的一種生成正態分布隨機數的方法進行分享-BOX MULLER法
BOX-MULLER 的Python實現
#import import numpy as np def boxmuller(): summa = 1 size = 1 x = np.random.uniform(size=size) y = np.random.uniform(size=size) z = np.sqrt(-2 * np.log(x)) * np.cos(2 * np.pi * y) q = z * summa return q print(boxmuller()[0])
本文對於BOX-MULLER不作詳細的介紹,僅作為拓展
拓展
隨機數好壞的判定
- 1.德聯安全辦的官方網址:
[BSI-STARTSEITE](傳送門)
- 2.NIST測試方法
[NIST.gov - Computer Security Division - Computer Security Resource Center](傳送門)
- 3.Diehard測試程序
[http://www.stat.fsu.edu/pub/diehard/](傳送門)
- 4.TestU01測試程序
[http://simul.iro.umontreal.ca/testu01/tu01.html](傳送門)
線性同余法
- [一種基於線性同余算法的偽隨機數產生器]論文 發表於[純粹數學與應用數學]21卷第三期(請自行尋找論文)
梅森旋轉法
- [線性反饋移位寄存器與梅森旋轉算法 - ACdreamer](傳送門)
- [wikipedia](傳送門)
Python下探究隨機數的產生方法