M06: integration

請務必詳閱作業描述 (一), (二), (三) 及 (四)

主講人: jserv / 課程討論區: 2024 年系統軟體課程
返回「Linux 核心設計/實作」課程進度表

Shannon Entropy 與亂度

「亂數」指結果不可預測的事件或數值序列。在統計學領域內，隨機事件的結果和在特定範圍內出現的機率，可藉由機率分布來展現。同時，隨機性也可透過事件發生時蘊含的資訊量進行衡量。Entropy (熵) 的主體概念是「資訊價值」，最初由 Claude Shannon 提出，他描述的資料傳輸系統由三大要素構成：來源資料、傳輸通道和接收端。通訊的根本問題在於探討接收端是否能夠僅憑傳輸通道輸出的資料準確推斷出來源資料。熵所代表的，是來源資料在經由傳輸通道壓縮時不遺失任何資訊的機率上限。

Information content 又稱為 self-information，反映我們在接收到某一訊息時獲得的資訊量。當某事件發生的機率達到 100% 時，該事件的資訊量將為零；反之，事件發生的機率愈低，其資訊量則愈高。因此，隨機事件

$x$ 的資訊本質僅與事件本身的發生機率相關。由於隨機事件是獨立的，且資訊本質僅與機率相關，資訊本質因此具備可加性。

觀察資訊本質，我們可歸納出兩個關鍵特性：

• 事件發生機率愈低，從資訊本質獲得的資訊量就愈高。
• 資訊本質具有可加性。

一個事件

$E$ 的 Information content 可定義為一個隨著

$E$ 發生機率

$p(E)$ 而下降的函數。由於對數函數

$\log$ 完美契合上述特性 (也僅有

$\log$ 函數能完美符合這些定義)，我們因此可以將資訊本質的函數定義為：

$$lo{g}_2(\frac{1}{p(E)})$$

其中，

$I$ 為資訊本體函數、

$P$ 為機率函數、

$x$ 為隨機事件，此定義符合

• 資訊本體與機率成反比，且當機率為 1 時，資訊本體為 0
• 兩隨機事件 a 與 b 的交集 c ，其資訊本體為
  $I(P(c))=I(P(a))+I(P(b))$

於是我們定義 information content 的公式為

$$I(E)=-lg(p(E))=lg(\frac{1}{p(E)})$$

以 Shannon entropy 來說，它代表需要問幾次是非題來猜中一個隨機值，定義如下

$$H(X)=-\sum _{\eta \in \chi }{p}^X(\eta )\cdot lg(p^X(\eta ))$$

$-lg(p^X(\eta ))$ 代表的即是字元

$\eta$ 的 informational content ，所以熵其實就是資訊本體的期望值，代表所接收的每則消息平均的資訊，亦可理解為亂度，當一則訊息越混亂，期熵越高。因此，可解讀如下

$$H(X)=-E[lg(p^X(\eta ))]$$

我們可由熵來計算亂度，由上面公式可知，熵最大的情況為，所有隨機事件都是平均分佈，也就是 uniform distribution ，此時

$H(X)=1/n\times P(x_1)+1/n\times P(x_2)+...+1/n\times P(x_n)$ 就是最大值。

〈由大數法則理解資訊熵的數學意義〉從資料壓縮編碼的視角探討 Shannon Entropy，這也正是 Shannon 最初提出他的理論時所採用的模型。透過大數法則的推導，找到一種方法來識別

$ϵ-$ high-probability set 的字串集合，發現要找到壓縮後錯誤概率小於

$ϵ$ 的字串集合，這樣的集合的總位元數約為

$2^{n\cdot H(S)}$，其中原始的每個字串長度為

$n$ 位元，因此壓縮比例為

$H(S)$。

Shannon Entropy 達到最大值的情形出現在

$P(s_1)=P(s_2)=\cdots =P(s_n)=\frac{1}{n}$ 時，即為均勻分佈（uniform distribution）。由於我們難以根據先前的資訊預測下一個出現的字元，這樣的字串更有可能被視為是隨機產生的。

熵值越大意味著包含的資訊量越豐富。在不損失原有資訊量的前提下進行壓縮，就必須使用更多位元。相反地，若熵值很小，則表示該字串集合易於預測，可以用更少的位元來表達原有的資訊，從而使壓縮帶來更大的效益。

英文字串的壓縮是個很好的案例，因為我們通常可根據最初幾個字元的組合預測後續將出現的字元，因此對英文字串進行壓縮效果顯著。在 Wikipedia Entropy 條目有個案例:

consider the transmission of sequences comprising the 4 characters 'A', 'B', 'C', and 'D' over a binary channel. If all 4 letters are equally likely (25%), one can not do better than using two bits to encode each letter. 'A' might code as '00', 'B' as '01', 'C' as '10', and 'D' as '11'. However, if the probabilities of each letter are unequal, say 'A' occurs with 70% probability, 'B' with 26%, and 'C' and 'D' with 2% each, one could assign variable length codes. In this case, 'A' would be coded as '0', 'B' as '10', 'C' as '110', and 'D' as '111'. With this representation, 70% of the time only one bit needs to be sent, 26% of the time two bits, and only 4% of the time 3 bits. On average, fewer than 2 bits are required since the entropy is lower (owing to the high prevalence of 'A' followed by 'B' – together 96% of characters).

假設今天有 4 個字元需要編碼，分別為 'A', 'B', 'C', 'D'， 且出現的機率分別為

$70\mathrm{\%}$,

$26\mathrm{\%}$,

$2\mathrm{\%}$,

$2\mathrm{\%}$，當我們依照這個機率建構出一個 optimal code tree 時：

當我們根據 optimal code tree 進行編碼時，將左側分支定義為 '0'，右側分支定義為 '1'。在該編碼規則下，'A' 被編碼為 '0'，'B' 為 '10'，'C' 為 '110'，'D' 為 '111'。這種設計考量在於 'A' 出現的機率最高，為了提高解碼效率，理應讓 'A' 的編碼盡可能短。相對地，由於 'C' 和 'D' 的出現機率較低，它們的編碼自然較長。

從熵的觀點分析，由於 '0' 和 '1' 所能承載的資訊量是等同的，我們可透過編碼長度來推斷各個字元的熵值。由於 'C' 和 'D' 的出現機率較小，它們蘊含的資訊量較大，因此對應的編碼也較長。然而，當 'A', 'B', 'C', 'D' 出現的機率均為

$0.25$ 時，此時的 optimal code Tree 如下：

4 個字元的編碼長度均為 2，optimal code tree 搜尋成本為

$2\times 0.25+2\times 0.25+2\times 0.25+2\times 0.25=2$，與前一案例的搜尋成本 1.34 相比，本例的搜尋成本較高，換言之，以熵值的觀點，相較於上個例子，本例的熵值較高。

當資料的熵越高時，我們需要足夠多的位元來攜帶這些資訊，也就較難壓縮，反之，當一組資料非常有規律時，這組資料的熵就會很低，我們就可用這些規律進行編碼與資料壓縮。

線性回饋移位暫存器

Linear feedback shift register (LFSR) 指遞迴關係能以線性函數 (具疊加性、一次齊次性的函數) 表達、主要以一串位移暫存器組成的電子電路，或上述電路所代表的演算法，其應用包括 PRNG 場景。

LFSR 範例:

1. 初始狀態
   $s=[s_2,s_1,s_0]=[1,0,0]$
2. 反饋
   $p=[p_2,p_1,p_0]=[0,1,1]$
3. 線性運算 s[n+1] = s[n]*p[2] ^ s[n-1]*p[1] ^ s[n-2]*p[0]

Fibonacci LFSR (又稱 Maximum-length LFSR) 由一列相接的移位暫存器組成，邏輯電平隨時脈朝右方傳遞，並透過對事先選定的幾個暫存器值做抑或運算決定最左方的數值。其值構成的序列稱為 maximum-length sequence (MLS, M-sequence）。

以軟體實作相應的 PRNG 時，能定義 lfsr 函式，其輸入值為指向無號整數的指標，該整數的各個位元即對應到電字電路各個移位暫存器的邏輯電平。以最高位對應上述最左端的暫存器，並以位元右移模擬移位暫存器的行為。

出處: Wikipedia

上圖展示一個 Fibonacci 線性回饋移位暫存器，以第 11, 13, 14, 16 位元的值做 XOR 運算後指派到第 1 位元，而第 k 位元的值指派到第 k + 1 位元（

$k\ge 1$）。

以下程式碼利用 ((*up) >> (64 - k)) & 1 取得

$2^{64-k}$ 位的值以做 XOR 運算，接著利用 ((*up) >> 1) | (new_bit << 63) 得到將運算結果指派到

$2^{63}$ 位，其他位右移一位的值，並指派回 *up 。

/* Implementation of LFSR (linear feedback shift register)
 * on uint64_t using irreducible polynomial x^64 + x^61 + x^34 + x^9 + 1
 * (On 32 bit we could use x^32 + x^22 + x^2 + x^1 + 1)
 */
static void lfsr(uint64_t *up)
{
    uint64_t new_bit =
        ((*up) ^ ((*up) >> 3) ^ ((*up) >> 30) ^ ((*up) >> 55)) & 1u;
    /* don't have to map '+1' in polynomial */
    *up = ((*up) >> 1) | (new_bit << 63);
    /* shift *up by 1 to RIGHT and insert new_bit at "empty" position */
}

多項式的乘方跟 64 位元整數第幾位數正好相反。第 0 位是

$x^{64}$ ，第 1 位是

$x^{63}$ 故式子中 (*up) >> 3 是代表

$x^{64-3}$ 為

$x^{61}$ 以此類推。Linux 核心原始程式碼 lib/random32.c 內有說明。

以線性回饋移位暫存器的原理實作的 PRNG 又稱為 Tausworthe 產生器 (Tausworthe generator），和 Robert C. Tausworthe 於 1965 年在《Mathematics of Computation》期刊發表的論文〈Random Numbers Generated by Linear Recurrence Modulo Two〉有關。這個別稱在 Linux 核心專案中也有出現。

執行

$ git log --grep='\(LFSR\|[Ll]inear [Ff]eedback [Ss]hift [Rr]egister\|[Tt]ausworthe\)'

能列出 Linux 核心原始碼中，所有提及 LFSR 的 commits 。截至 v6.4-rc6 （commit 858fd168a95c），在 master 分支上共有 23 個 commits 含有此關鍵字。再以關鍵字搜尋程式碼，能幫助我們了解 LFSR 實作於 Liunx 核心的過程。

Linux 核心專案於 2.6.2-rc2 開始採用 Git 版本控制系統，對應的 commit 1474855d0878 也是專案中最早的 commit 。在這之前，Fibonacci LFSR 實作就已出現於 drivers/net/ewrk3.c 裡頭的 get_hw_addr 函式中，作為乙太網路卡裝置驅動程式的一部分。

另外在 net/core/utils.c 當中，net-random 函式實作 taus88 演算法，是一種 maximally equidistributed combined Tausworthe generator，是由 Pierre L'Ecuyer 設計的演算法，週期大約為

$2^{88}$。

當時雖然已有 include/linux/random.h ，但直到 commit aaa248f6c9c8 （發布於 2.6.19-rc3 ），net_random 才獨立成 lib/random32.c 當中的 random32，其原型宣告於 include/linux/random.h 中。

random32 經過一系列改進，包含 reseeding 與升級演算法到 taus113 （commit a98814cef879）。

Fibonacci LFSR 或 Galoise LFSR 仍被當做一種簡單的 PRNG ，在 Linux 核心的 driver 和 profiler 等子系統採納。。

亂數產生器

以下由 Shawn531 提供

前述 Entropy 的極大值發生在

$P(x_1)=P(x_2)=...=P(x_n)=\frac{1}{n}$ 時

$(S=lo{g}_{b}n)$， 此分佈即爲 uniform distribution。因此我們要設計一個亂數產生器，讓它趨近於 uniform distribution 的分佈。Pseudorandom_number (PRN) 是透過一些手段讓亂數雖然是以確定的方式產生，但卻可以在我們較小的使用範圍內不可預測，PRNG 即為 PRN 的產生器。

Linux 核心提供 /dev/random 作為 PRNG。將產生的亂數提供給 ent 工具。

$ head -c 10M /dev/random | ent
Entropy = 7.999982 bits per byte.

Optimum compression would reduce the size
of this 10485760 byte file by 0 percent.

Chi square distribution for 10485760 samples is 255.57, and randomly
would exceed this value 47.82 percent of the times.

Arithmetic mean value of data bytes is 127.4605 (127.5 = random).
Monte Carlo value for Pi is 3.144157846 (error 0.08 percent).
Serial correlation coefficient is 0.000057 (totally uncorrelated = 0.0).

以下比較xorshift64, xorshift128p, xoshiro256+, xoshiro256** 等 PRNG 的亂度。

xorshift64

Xorshift 是 Linear-feedback shift register (LFSR) 的子集合，簡單透過使用左移右移和 xor 對先前的狀態進行運算，xorshift 有若干變形，各種位元長度的變形。xorshift32 會有一個 32 位元的狀態，xorshift64 會有一個 64 位元的狀態，xorshift128會有 4 個 128 位元的狀態。要產生隨機數很重要的一點是初始化 state ，且初始化不能為 0。

	uint64_t x = state->a;
	x ^= x << 13;
	x ^= x >> 7;
	x ^= x << 17;
	return state->a = x;

xorshift128p

xorshift128p 會在最後透過加法達成非線性的轉換，相比在最後加上乘法的 xorshift128star。

	uint64_t t = state->x[0];
	uint64_t const s = state->x[1];
	state->x[0] = s;
	t ^= t << 23;		// a
	t ^= t >> 18;		// b -- Again, the shifts and the multipliers are tunable
	t ^= s ^ (s >> 5);	// c
	state->x[1] = t;
	return t + s;

xoshiro256p 及 xoshiro256ss

xoshiro256p 和 xoshiro256ss在架構相同，差別在在非線性項的產生。用於 GNU Fortran compiler 和 Lua 專案。

uint64_t rol64(uint64_t x, int k)
{
    return (x << k) | (x >> (64 - k));
}

uint64_t xoshiro256p(struct xoshiro256p_state *state)
{
    uint64_t *s = state->s;
    uint64_t const result = s[0] + s[3];
    uint64_t const t = s[1] << 17;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;
    s[3] = rol64(s[3], 45);

    return result;
}

初始化狀態時，建議用與初始產生器完全不同且不會得到全 0 的產生器，進行初始化。以下使用 slpitmix64 位的 seed 產生器。在實作部分需要針對每種 PRNG 做不同狀態的初始化。

uint64_t splitmix64(struct splitmix64_state *state)
{
    uint64_t result = (state->s += 0x9E3779B97f4A7C15);
    result = (result ^ (result >> 30)) * 0xBF58476D1CE4E5B9;
    result = (result ^ (result >> 27)) * 0x94D049BB133111EB;
    return result ^ (result >> 31);
}

若狀態為s[4]，我們需要對 4 位分別做 splitmix64 的動作。之後就可以直接呼叫產生器。

void xorshift128_seed(const uint64_t seed)
{
    seeded = true;
    splitmix64_seed(seed);
    for (int i = 0; i < 2; i++)
        s[i] = splitmix64();
}

commit 2331493

在整合的 qrandom.c 中，0 為系統內建的 PRNG，1 為 xoshiro+，2 為 xoshiro**，3 為xorshift64，4 為xorshift128+。會選用這幾種 PRNG 就是因為他們的 state 皆是 uint64_t 的資料型態，不同的就只有 state 個數，方便使用 splitmix64 實作初始化。

實驗結果

PRNG	內建	xorshift64	xorshift128+	xoshiro+	xoshiro**
Entropy	7.999982	7.999981	7.999982	7.999983	7.999984
Compression	0%	0%	0%	0%	0%
Chi square distribution	258.69	269.09	255.15	244.50	227.94
confidence level	42.39%	26.04%	48.56%	67.08	88.75%
Arithmetic mean vlue	127.4907	127.4756	127.5119	127.5036	127.4449
Monte Carlo value	3.143301828	3.142072732	3.141436440	3.142024667	3.143034036
Serial correlation coefficient	-0.000258	0.000169	0.000142	-0.000350	0.000024

1 個位元組的字元能夠表示的空間為

$2^8$，則最大的 entropy 為

$S={\log }_2(2^8)=8$，可見這次測試的 5 種 PRNG 的 entropy 皆非常接近最大值 8，最大壓縮率皆趨近於 0。算術平均、蒙地卡羅法計算 Pi 值、及序列相關係數在 5 種 PRNG 皆相當接近，而對照參考值顯示五種 PRNG 皆足夠分散。

這一次實驗最能看出 5 種 PRNG 的效能是 Chi square distribution，信心水準離 5% 越遠，代表產生的亂數越可能符合 uniform distribution，xorshift64 為最差，這是由於只有他是線性的，其他的 PRNG 都是非線性的變形。xoshiro** 為這一次表現最好的 PRNG，信心水準來到了 88.75%，可能會是常態分布，但由於 xoshiro** 的運算時間會比 xoshiro+ 來的多，因此權衡之際下xoshiro+ 已成為相當多框架的 PRNG。

xoroshiro 亂數產生器

xoroshiro 的命名由 xor, rotate, shift, rotate 等詞彙組合，是 Xorshift 的改良。尋常的位移操作從一端溢出的位元會被丟棄，但 rotate 操作會讓它們從另一端出現。
若對 1011 進行向左移 2 個位元操作，會得到 1100，但若是 rotate 操作，則會得到 1110。

xoroshiro128+ 擁有 128 位元的內部狀態，分為 2 個 64 位元，並返回一個 64 位元的隨機數值。若內部狀態由 64 位元整數

$X$,

$Y$ 組成，則透過

$X+Y$ 獲取隨機數。

狀態轉移藉由以下過程進行：

X = rotl(X, 24) ^ (X^Y) ^ ((X^Y) << 16)
Y = rotl(X^Y, 37)

計算過程：
rotl(X,n) =:

$R_l(n)\in M_{64}({\mathbf{F}}_2),$

$$\begin{array}{rl}\left[\begin{array}{c}{x}_{n+1}\\ y_{n+1}\end{array}\right]& =\left[\begin{array}{c}{R}_l(24)x_n+x_n+y_n+L_{16}(x_n+y_n)\\ R_l(37)(x_y+y_n)\end{array}\right]\\ & =\left[\begin{array}{c}(R_l(24)+E+L_{16})x_n+(E+L_{16})y_n\\ R_l(37)x_y+R_l(37)y_n\end{array}\right]\\ & =\left[\begin{array}{cc}{R}_l(24)+E+L_{16}& E+L_{16}\\ R_l(37)& R_l(37)\end{array}\right]\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right].\end{array}$$

$R_l(n)\in G{L}_{64}({\mathbf{F}}_2)$ 於是

$R_l(n{)}^{-1}=R_l(64-n)$，對

$L_{16}^{\prime }:=E+L_{16}$ 也是如此。在

${\mathbf{F}}_2^{64}$ 中，

$b=a+(a<<n)⟺a=b+(a<<n)$，當

$x,y$ 的第

$i$ 個位元分別為

$x_i,y_i$ 時，

$$x_i=\{\begin{array}{l}{y}_i(1\le i<n)\\ y_i+x_{i-n}(n\le i\le 64)\end{array}$$

由於它是可逆的，亦即我們得到逆矩陣。為了簡單起見，假設

$R_l(24)=A,L_{16}^{\prime }=B,R_l(37)=C$，於是：

$$\left[\begin{array}{c}{x}_{n+1}\\ y_{n+1}\end{array}\right]=\left[\begin{array}{cc}A+B& B\\ C& C\end{array}\right]\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]$$

現在求其逆矩陣：

$$\left[\begin{array}{cc}A+B& B\\ C& C\end{array}\right]\left[\begin{array}{cc}{X}_1& X_2\\ X_3& X_4\end{array}\right]=\left[\begin{array}{cc}E& 0\\ 0& E\end{array}\right]$$

於是：

$\begin{array}{rl}(A+B)X_1+B{X}_3& =E\\ (A+B)X_2+B{X}_4& =0\\ C{X}_1+C{X}_3& =0\\ C{X}_2+C{X}_4& =E\end{array}$

上述 4 式成立。

$A,B,C\in M_{64}({\mathbf{F}}_2)$ 及

$-X=X$，注意到乘法是不可交換的，基於這點求解。根據下面的 2 個式子：

$\begin{array}{r}{X}_3=X_1\\ X_4=C^{-1}+X_1.\end{array}$

將此依次代入到上面的 2 個式子中。首先從第一個開始:

$\begin{array}{rl}(A+B)X_1+B{X}_1& =E\\ A{X}_1+B{X}_1+B{X}_1& =E\\ A{X}_1& =E\\ X_1& =A^{-1}.\end{array}$

接著是：

$\begin{array}{rl}(A+B)X_2+B(C^{-1}+X_2)& =0\\ A{X}_2+B{X}_2+B{C}^{-1}+B{X}_2& =0\\ A{X}_2+B{C}^{-1}& =0\\ A{X}_2& =B{C}^{-1}\\ X_2& =A^{-1}B{C}^{-1}.\end{array}$

從上可知：

$\begin{array}{rl}{\left[\begin{array}{cc}A+B& B\\ C& C\end{array}\right]}^{-1}& =\left[\begin{array}{cc}{A}^{-1}& A^{-1}B{C}^{-1}\\ A^{-1}& A^{-1}B{C}^{-1}+C^{-1}\end{array}\right]\\ & =\left[\begin{array}{cc}{R}_l(24{)}^{-1}& R_l(24{)}^{-1}{L}_{16}^{\prime }{R}_l(37{)}^{-1}\\ R_l(24{)}^{-1}& R_l(24{)}^{-1}{L}_{16}^{\prime }{R}_l(37{)}^{-1}+R_l(37{)}^{-1}\end{array}\right]\\ & =\left[\begin{array}{cc}{R}_l(40)& R_l(40)L_{16}^{\prime }{R}_l(27)\\ R_l(40)& R_l(40)L_{16}^{\prime }{R}_l(27)+R_l(27)\end{array}\right].\end{array}$

不需要

$L_{16}^{\prime }$ 的可逆性。

也就是說：

$$\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]=\left[\begin{array}{cc}{R}_l(40)& R_l(40)L_{16}^{\prime }{R}_l(27)\\ R_l(40)& R_l(40)L_{16}^{\prime }{R}_l(27)+R_l(27)\end{array}\right]\left[\begin{array}{c}{x}_{n+1}\\ y_{n+1}\end{array}\right]$$

$$\begin{array}{rl}{x}_n& =R_l(40)x_{n+1}+R_l(40)(E+L_{16})R_l(27)y_{n+1}\\ & =R_l(40)x_{n+1}+R_l(3)y_{n+1}+R_l(40)L_{16}{R}_l(27)y_{n+1}\\ y_n& =R_l(40)x_{n+1}+R_l(40)(E+L_{16})R_l(27)y_{n+1}+R_l(27)y_{n+1}\\ & =R_l(40)x_{n+1}+R_l(3)y_{n+1}+R_l(40)L_{16}{R}_l(27)y_{n+1}+R_l(27)y_{n+1}.\end{array}$$

綜合以上，可知：

X = rotl(X,40) ^ rotl(Y,3) ^ rotl((rotl(Y,27)<<16),40)
Y = rotl(X,40) ^ rotl(Y,3) ^ rotl((rotl(Y,27)<<16),40) ^ rotl(Y,27)

回顧作業一的洗牌過程

以下由 weihsinyeh 貢獻

《The Art of Computer Programming》提及 shuffle 過程的四步驟 :

1. Initialize : 是將序列的數量存給 j 。
2. Generate
   $U$ : 產生一個介於 0~1 的數字 。
3. Exchange : Set
   $k$ <-
   $\lfloor j\times U\rfloor +1$. Exchange
   $X_j$ <->
   $X_k$ 。
4. Decrease
   $j$ : if
   $j>1$ 回到第二步。

每次沒被選過的都是 1~j，因為每次都會選一個節點交換，而選的那個節點不需要是隨機的，只要交換一個節點為隨機的便可以達到 shuffle 的功能，因此才將另一個節點都是

$j$。因此我將演算法修正用迴圈每回合減少 j，且每回合選一個值跟

$j$ 交換。

comit 88ae9fd

至於取隨機數時，究竟該產生一個介於 0 到 1 的數值再做乘法，還是用取餘數的方式呢？隨機數的範圍都只有 0 到 RAND_MAX。在 C99 標準中說明 "RAND_MAX assumed to be 32767"，這代表隨機數的範圍會被限制。

int rand(void)  // RAND_MAX assumed to be 32767
{
    next = next * 1103515245 + 12345;
    return (unsigned int)(next/65536) % 32768;
}

C99/C11 Standard (§ 7.23.2.4) The time function - Returns
The time function returns the implementation's best approximation to the current calendar time. The value (time_t)(-1) is returned if the calendar time is not available. If timer is not a null pointer, the return value is also assigned to the object it points to.

srand(time(NULL)) 將目前的時間的放進 next。傳遞參數 NULL 代表沒有要將值存到一個指標。但這沒有解決隨機數的範圍會被限制的問題。同時，終於找到如何取介於 0 到 1 的浮點數值。首先發現這裡對 RAND_MAX 取餘數，所以只要再除 RAND_MAX 的大小就可以得到 0 到 1 間的浮點數值。

而

$2^{16}=65536$ 、

$2^{15}=32768$ 這二個數字一定跟 next 的型態有關連，且unsigned long int 的 minimum size (bits) 是 32。因為

$0$~

$2^{32}$ 除

$2^{16}$ 得到的範圍在

$0$~

$2^{16}$ 之間。範圍只取　

$2^{15}$ 則是因為可以轉型到 unsigned int。也就是說，若要得到 0 到 1 浮點數值，可用 double result = rand()/(RAND_MAX+0.0);。

    srand(time(NULL));
    for (; len != 0 ;safe = safe->prev) {
        int random = rand() % (len--);

srand(time(NULL)) 只初始化以取得當下時間一次，但每次卻都得到隨機的數字。其實初始化 next 這個變數後，就不再跟時間有關，不會再去呼叫 time() 拿取當下的時間。會動的只有 next ，而變化的關係是乘 1103515245 再加 12345。因此，這個亂數有規律，一點都不亂。所謂的 random seed 一點也不是產生亂數的由來，這個名字誤導我很久。它其實是找一個獨一無二的數字，也就是當下的時間作為一個起點。讓下次若要再 srand(time(NULL)) ，會產生像隨機的序例，因為他們的起點不一樣。

The same initializing value of seed will always return the same random sequence.

接下來是乘 1103515245 與加 12345 的由來：
這源自於 Linear congruential generator (LCG)。

$X_{n+1}=(a{X}_n+c)modm$。從上面的例子明顯看得出來這有規律，會發生重複，即便在最好的情況也會在

$m$ 次後就出現重複，而重複的週期取決於

$m$、

$a$、

$c$ 的值，一個數值改變，週期就會變。因此為了產生看起來像是亂數的序列，就要去選擇適當的

$m$、

$a$、

$c$ 讓重複的頻率將低，也就是要增加週期的長度 period length 。

參考《Numerical Recipes in C》的解釋

$c$ 是

$2^{15}$ 的原因：

arithmetic done on unsigned long quantities is guaranteed return the correct low-order bits

選擇的理論來自《The Art of Computer Programming》3.2.1 The linear congruential sequence

$m$、

$a$、

$c$ and

$X_0$ has period length m if and only if

1. c is relatively prime to m;
2. b = a − 1 is a multiple of p, for every prime p dividing m;
3. b is a multiple of 4, if m is a multiple of 4.

統計方法驗證 shuffle

樣本數量為 4 ，自由度為 4!= 24 - 1 。

p-value (probability)	0.995	0.99	0.975	0.95	0.90	0.10	0.05	0.025	0.01	0.005
freedom : 23	9.260	10.196	11.689	13.091	14.848	32.007	35.172	38.076	41.638	44.181

• 實驗ㄧ

用統計方法的驗證，很明顯完全不符合 Uniform Distribution 。

chi square sum: 36487.09585753372 其對應的 p value < 0.05 。不支持 uniform distribution 。

• 實驗二

修正程式碼用統計方法的驗證，很明顯完全跟作業說明一樣符合 Uniform Distribution。

修正方式是把原先在上方圖中的 srand(time(NULL)); 這行移掉。它的功能是每次 shuffle 都會將 Linear congruential generator 的 seed 初始化為當下的時間。但這樣其實每次shuffle 都是一樣的操作，根本沒有shuffle 到，因為每次 seed 都被 static unsigned long int next = 1; 這樣初始化。為其實 shuffle 是有規律的，只是不知道它從週期的哪個起點開始。那有規律不就不是獨立事件。

chi square sum: 24.600921614745832 其對應的 p value > 0.05 。支持 uniform distribution 。

但再做一次實驗發現他們的數值又不一樣同時也沒有規律，但每次都確實是用 next = 1，重新開始。所以我又寫一個程式去測試最基本的東西，我發現不加初始化的時候，用一個迴圈去跑 rand 。每次重新執行程式的確跑出來的是會一樣，這與上面的理解一樣。

所以代表 shuffle 的程式碼有一個會變化的東西。我發現原來是在 q_test.c 裡面有這 srand(os_random(getpid() ^ getppid())); 所以才讓每次都會初始化。但用這種方式 os_random 初始化，跟用時間去初始化竟然效果差如此多。

但要確認兩件事，到底是每次都在 q_shuffle 裡面做 srand(time(NULL)); 才導致上面的第一張圖沒有 Uniform Distribution 。還是因為改成 srand(os_random(getpid() ^ getppid())); 不用時間初始化的原因。

但由於我認為程式執行一次期間每次 q_shuffle 的這件事情應該需要跟上一次沒有關係。若每次都沒有在 q_shuffle 初始化 seed ，其實從程式執行一開始，就可以知道 shuffle 完的所有結果，因為我只有一開始去變化 seed ，後面一切操作都是有跡可尋的，就只是按照一定的規則在交換佇列裡面元素的位置而已。

• 實驗三

將 q_test.c 裡面的初始化改為 srand(time(NULL)); 。

chi square sum: 12.213699419190709 其對應的 p value > 0.05 。支持 uniform distribution 。
在 q_test.c 的 main 函式裡初始化只是讓每次程式重跑一次會得到不一樣的結果而已。 所以用 srand(os_random(getpid() ^ getppid())); 還是 srand(time(NULL)); 都可以得到 Uniform Distribution 。

• 實驗四

在每次呼叫 q_shuffle 函式前先有初始化 seed 的步驟，實作方式是在 do_shuffle 裡面增加 srand(os_random(getpid() ^ getppid()));

chi square sum: 18.36418182690923 其對應的 p value > 0.05 。支持 uniform distribution 。

原先還沒做實驗前我以為是 q_shuffle 裡面「有無」加入初始化 seed 的步驟影響很大。但做完實驗才發現是srand(os_random(getpid() ^ getppid())); 與srand(time(NULL));差別。因此我先前是不需要去懷疑這個實作方式能夠說明是 uniform distribution 的，是我自己沒有用對的方式實作，而第二個作法即便最後呈現的是 uniform distribution 的樣子，但它仍然是錯的，因為不是 q_suffle實作才導致呈現成 uniform distribution ，而是 Linear congruential generator 本身在週期內不重複的特性才得以使第二個作法呈現 uniform distribution 。此外由於每一次 rand() 都是的 next 都跟前一次 next 有關係，所以每次 shuffle 都不是獨立事件，這次的 shuffle 影響下一次 shuffle ，因此根本不符合 uniform distribution 。

第四個作法就解決此問題，每次 shuffle 前都透過srand(os_random(getpid() ^ getppid())); 初始化。讓 shuffle 不影響下次 shuffle 為獨立事件，符合 uniform distribution 。所以重點是 random seed 的作法。

接下來看 commit f17a9a8 變動 seed 的方法。步驟一、先用 os_random 函式的位址對傳進去的　seed = (getpid() ^ getppid()) 作 xor 。而這個操作完，不只用每一次執行一次程式就透過 pid 與 ppid 得到，函式的位址也是 Address space layout randomization 隨機化過得。但想了一下 os_random 的函式其中接下來的步驟二、是取得目前的高解析度時間 (high-resolution time)，然後步驟三、再做 random_shuffle 。因為 ASLR 在程式碼一執行後函式的位址就不再變動，所以他的功能是防止 exploitation of memory corruption vulnerabilities。

那讓 shuffle 變成獨立事件的就是後面步驟二與步驟三。所以我就先不做步驟一看結果會不會影響到亂數這件事情，也就是把程式碼的 x = seed 直接做後面步驟。實驗完發現不影響仍然是 uniform distribution，所以步驟一的確主要的功能是避免漏洞。

開發 random_shuffle 的時候，參考〈Prospecting for Hash Functions〉，讓操作不可逆。作者 Chris Wellons 發現達成不可逆的兩個性質，他尋找適當的雜湊函數，藉由程式 prospect 一直去算看是否有盡可能滿足兩個特性 High avalanche effect 與 Low bias。 而 prospect 自己算出一個很好的 pattern :

The pattern is XOR-right-shift, multiply, XOR-right-shift, multiply, XOR-right-shift

這個操作因是常數時間，不常被使用在密碼學。

同時用這三個運算子後的結果仍需是 1-to-1 mapping 的函數，否則又會有分佈不均勻的問題。

The XOR-right-shift diffuses bits rightward and the multiply diffuses bits leftward. I like to think it’s “sloshing” the bits right, left, right, left.

該操作的概念：比如原先是 ABCD

$\oplus$ CDAB ，操作後得到的資訊是（A 跟 C）與（B 跟 D）兩組是同或異。XOR 的結果會導致 higher bytes 與 lower bytes 是一樣的，這對隨機這件事情沒有幫助。接下去看 multiply xor shift 裡面解釋 shift 將高位元移到低位元的原因：高位元會因低位元乘法進位導致變化大。而這個操作要做至少兩次，因為要讓 shift 到低位元「原先的高位元」也影響高位元「原先的低位元」（再乘法進位一次）。如此一來才能達成讓高低位元變化量平衡的目的。

之所以要做步驟三是因為又將取得的當下時間變成一個

It's statistically indistinguishable from a random permutation of all 32-bit integers.

而這個目的是為了讓隨機數不可逆。

接下來把步驟三拿掉看結果。發現移掉也是 uniform distribution 。所以第三步驟也不是讓每次執行都變獨立事件的關鍵，於是聚焦在第二步驟。所以現在就要實驗第二步驟改成 time(NULL) 其他不變的結果。做完實驗後發現真的是第二步驟讓整建事情變成獨立事件。

原來取得目前時間的方法會讓結果差很大。由 high-resolution time) 得知原先用 time(NULL) 初始化 seed 是用 wall clock's unsafe current time，改成 monotonic clock's unsafe current time。一如命名所及，monotonic clock 只存在一個 process 可以看到。這讓我想起以前第一次看到 tv_secv 是在寫 context switch 的時候，所以不像 wall clock 是可以被多個行程共用。也因此獨立事件這個說法有問題，使用 monotonic clock 是讓取時間這個操作變得更均勻分佈，因爲每次 shuffle 的因為迴圈會跑 sample 數量的次數，每輪 rand 再交換，操作所花費為 constant time。所以第二步驟中每次會依時間變化的變數更均勻的分佈。

attacker may be able to obtain high-resolution estimates through repeat execution and statistical analysis.

然而 monotonic clock 得到更精準的時間，讓變數變化均勻分佈。原先我都認為要得到亂數，不是越不規律越好，像前面提到的 LCG 就是因為存在規律，所以才要用 random seed。摘錄 Prospecting for Hash Functions：

One of the important properties of an integer hash function is that it maps its inputs to outputs 1:1. In other words, there are no collisions. If there’s a collision, then some outputs aren't possible, and the function isn’t making efficient use of its entropy.

雖然越不規律較好，但不規律的前提就是要變數均勻分布，否則若一個樣本出現次數過高，那就不能叫做亂數。所以一開始才用 Pearson's chi-squared test 檢驗虛無假說（Null hypothesis），用這個統計方法驗證是否 shuffle，此正好解釋為什麼用 monotonic clock 才能產生不拒絕虛無假設的 uniform distribution。然而因為時間精準度提高的同時又有 timing side-channel attacks 的問題。因此從此也知曉做第三步驟的原因，為了避免攻擊者透過推算出時間做攻擊，因此要讓第二步驟的變數不可逆推。

這整件事情來看每步驟的功能。第一步驟是為了避免漏洞，所以搭配 ASLR。第二步驟為了讓每次操作的變數符合虛無假說，所有樣本出現的機率相同使用高精度。第三步驟的意義在於，因為一些系統上面 ASLR 也無法杜絕漏洞被攻擊，同時還有第二步驟會導致 side-channel attack 發生，所以還要第三步驟的 random-suffle 使變數不可被攻擊者逆推。 commit 7a51365

延伸閱讀：Uniting the Linux random-number devices，亂數產生器的品質是近期 Linux 核心的重要開發方針之一。

xoro 核心模組

ksort 提供一個名為 xoro 的 Linux 核心模組，提供 /dev/xoro 裝置檔案，讓使用者得以取得 xoroshiro128+ 產生的亂數，參見 test_xoro.c 的使用。