L01: lab0

主講人: jserv / 課程討論區: 2023 年系統軟體課程

Image Not Showing Possible Reasons
The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported
Learn More →

返回「Linux 核心設計/實作」課程進度表

在 qtest 提供新的命令 `shuffle`

利用 Fisher–Yates shuffle 演算法來實作洗牌（shuffle）

先用 q_size 取得 queue 的大小 len。
隨機從 0 ~ (len - 1) 中抽出一個數字 random，然後 old 將指向從前面數來第 random 個節點，new 會指向最後一個未被抽到的節點，將 old 和 new 指向的節點的值交換，再將 len - 1。
隨著 len 大小變小，已經被抽到過，並交換值到 queue 後面的會愈來愈多，直到所有的節點都已經被抽到過，shuffle 就結束。

統計方法驗證 `shuffle`

Pearson's chi-squared test 能檢驗虛無假說（Null hypothesis），即某特定事件在樣本中觀察到的頻率分佈與特定理論分佈一致，事件必須是互斥且機率為 1。

如果要 shuffle 三個不同的數字，會出現六種可能的結果，彼此是互斥的，且發生機率加起來為 1。那假設 shuffle 是公平的（六種結果發生的機率相同），並遵守 Uniform distribution，那虛無假說為：

$H_{0}$ （虛無假說）: shuffle 的結果發生的機率相同，遵守 Uniform distribution
$H_{1}$ （對立假說）: shuffle 的結果發生的機率至少有一個不同

接下來透過 Pearson's chi-squared test 來驗證 shuffle 三個數字的頻率是符合 Uniform distribution:

1. 先計算 chi-squared test statistic
$X^{2}$

X^{2} = \sum_{i = 1}^{n} \frac{(O_{i} - E_{i})^{2}}{E_{i}}

$X^{2}$ : Pearson's cumulative test statistic
$O_{i}$ : the number of observations of type i
$E_{i}$ : the expected count of type i

E:  166666
O:  {'abc': 166391, 'bac': 166829, 'bca': 166807, 'acb': 166790, 'cab': 166862, 'cba': 166321}
0.45375181500726003
0.15941463765855063
0.11928647714590858
0.0922563690254761
0.23049692198768795
0.7141528566114265
sum:  1.76935907743631

在測試 shuffle 1000000 次後，六種結果各自的頻率如下表：

	觀察到的頻率	期待的頻率	$\frac{(O_{i} - E_{i})^{2}}{E_{i}}$
[1, 2, 3]	166391	166666	0.45375181500726003
[2, 1, 3]	166829	166666	0.15941463765855063
[2, 3, 1]	166807	166666	0.11928647714590858
[1, 3, 2]	166790	166666	0.0922563690254761
[3, 1, 2]	166862	166666	0.23049692198768795
[3, 2, 1]	166321	166666	0.7141528566114265
Total			1.76935907743631

X^{2}

= 1.76935907743631

2. 決定自由度（degrees of freedom）

對於 N 個隨機樣本而言，自由度為 N - 1。我們 shuffle 3 個數字會有六種結果，因為加起來機率為 1 ，所以實際上可以自由變換的變數只有五個，其中一個結果的機率為 1 減去另外五個結果發生的機率，所以自由度為 5。

3. 選擇顯著水準

顯著水準（Significance level）α 代表在虛無假說（
$H_{0}$ ）為真下，錯誤地拒絕
$H_{0}$ 的機率，即型一錯誤發生之機率。
α = P(拒絕
$H_{0}$ |
$H_{0}$ 為真)
我們 alpha 設定為常見的 0.05。
P value
從卡方分布表找合適的 P value，我們的自由度為 5，
$X^{2}$ 為 1.76935907743631，因為 1.61 < 1.76935907743631 < 2.34，於是 P value 介於 0.9 和 0.8 之間。
Image Not Showing Possible Reasons
- The image was uploaded to a note which you don't have access to
- The note which the image was originally uploaded to has been deleted
Learn More →

4. 統計檢定的結果不拒絕虛無假說

對於要拒絕的虛無假說（

H_{0}

），觀察到的結果必須具有統計顯著性，即觀察到的 P value 小於預先指定的顯著水準 α。

因為 P value（0.8~0.9）> alpha（0.05），統計檢定的結果不拒絕虛無假說（

H_{0}

）
也就是樣本資料不拒絕 shuffle 的結果遵循 Uniform distribution，因為沒有足夠證據推翻。

從圖表來看 shuffle 的頻率是大致符合 Uniform distribution 的。

Image Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

測試程式

scripts/driver.py 用到 subprocess.call() 來產生新行程以執行指定的命令，然後等待命令完成後取得 return code，後續再檢查 return code 是否為 0 來得知測試是否成功。

因為這裡的測試我想取得 stdout，所以改用 subprocess.run() ，會執行參數中指定的命令，然後等待命令完成後，會回傳一個 CompletedProcess instance。
取得 CompletedProcess.stdout 再做編碼後，可以得到 stdout 的字串。
後續再將得到的輸出字串取出 shuffle 過後的結果，放到 nums 變數中。

測試用的 Python 腳本

import subprocess
import re
import random
from itertools import permutations
import random

# 測試 shuffle 次數
test_count = 1000000
input = "new\nit 1\nit 2\nit 3\nit 4\n"
for i in range(test_count):
    input += "shuffle\n"
input += "free\nquit\n"

# 取得 stdout 的 shuffle 結果
command='./qtest -v 3'
clist = command.split()
completedProcess = subprocess.run(clist, capture_output=True, text=True, input=input)
s = completedProcess.stdout
startIdx = s.find("l = [1 2 3 4]") 
endIdx = s.find("l = NULL")
s = s[startIdx + 14 : endIdx]
Regex = re.compile(r'\d \d \d \d')
result = Regex.findall(s)

def permute(nums):
    nums=list(permutations(nums,len(nums)))
    return nums

def chiSquared(observation, expectation):
    return ((observation - expectation) ** 2) / expectation 

# shuffle 的所有結果   
nums = []
for i in result:
    nums.append(i.split())

# 找出全部的排序可能
counterSet = {}
shuffle_array = ['1', '2', '3', '4']
s = permute(shuffle_array)

# 初始化 counterSet
for i in range(len(s)):
    w = ''.join(s[i])
    counterSet[w] = 0

# 計算每一種 shuffle 結果的數量
for num in nums:
    permutation = ''.join(num)
    counterSet[permutation] += 1
        
# 計算 chiSquare sum
expectation = test_count // len(s)
c = counterSet.values()
chiSquaredSum = 0
for i in c:
    chiSquaredSum += chiSquared(i, expectation)
print("Expectation: ", expectation)
print("Observation: ", counterSet)
print("chi square sum: ", chiSquaredSum)

產生直方圖的 Python 程式：

import matplotlib.pyplot as plt
import numpy as np
permutations = counterSet.keys()
counts = counterSet.values()

x = np.arange(len(counts))
plt.bar(x, counts)
plt.xticks(x, permutations)
plt.xlabel('permutations')
plt.ylabel('counts')
plt.title('Shuffle result')
plt.show()

測試 [1, 2, 3, 4] 做 shuffle 1,000,000 次的直方圖：

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

亂數

討論 Random Number Generator 夠不夠亂之前，要先回答以下二個問題：

What is random?
How to measure the randomness?

亂數直覺而言是個很簡單的概念，一般來說，我們認爲一個事件「亂」，表示我們無法預測該事件的結果，但統計學卻讓所謂的「亂」有跡可循，透過機率分佈的概念，我們即便無法準確的預測下一次的事件結果爲何，卻可以很「確定」的表示隨機事件的結果落在某個範圍內的機率是多少，或長期而言結果的統計量會趨近於某個確定的值 – 我們可以用機率分佈來「描述」隨機事件。

另一個角度來看，我們可能認爲一段文字是隨機的，如果我們無法從該段文字中讀出資訊，例如 "woif ejxv nwe" 可能是一個隨機的字串而 "what the f*ck" 則不是隨機 – 資訊的含量也可以用來描述隨機性。

從上面的兩個觀點中，我們可以截取出兩個概念，並討論看看這他們與亂數的關係 – probability & information

一個事件的結果無跡可循，換句話說結果發生在任何位置的機率是均等的，這樣的性質可以用 uniform distribution 來表達，但光是 uniform 是不夠的，考慮一個數列:

f (i) = a + (i m o d (b - a)), i = 1, 2, 3, . . .

可以預期結果會是分佈在

[a, b]

的正整數，且發生的頻率均等，但我們不會認爲這樣的數列是隨機的，因爲我們可以很容易的發現事件之間有前後關係，所以還要再多加上一個條件，事件發生的機率互相獨立。

在此觀點下，理想的亂數應符合以下分佈

X \sim^{i i d} U [a, b]

那麼以文字的資訊量切入的話什麼是理想的亂數呢？

我們得先談論何爲「資訊」(information)，考慮剛才的例子，爲何我們有辦法解讀出 "what the f*ck" 的含義？首先從單字來看，我們知道 what 按順序出現時可表達出某個意思，換句話說當出現 wha 時，有很高的機率接下來會出現 t；從句子來看，當我們看到 what 知道他是疑問詞，根據文法結構可以預期接下來會出現名詞的機率很高，相反如果出現的不是名詞我們將難以判斷該文句的意義，我們就認爲該句子是隨機而無意義的。

上述過程我們可以發現一件事情，當我們預期某個事件發生的機率很高，而他真的發生的時候，我們獲得的資訊其實沒有變多，也就是說上面這句話裏面事實上有許多部分是多餘的，即便我們將一部分的字元移除仍可以傳達相等的資訊（缺失的部分可以很容易預測），例如 "wat th fk" 即使沒有一個單字拼對我們仍可讀出一樣的訊息，兩個句子每個字元平均所含的資訊量是不相等的，這也是資訊壓縮的基本原理，要怎麼用最少的字元數表達相同的資訊量。

總結來說，當一個事件發生的機率越低，而他真的發生時，我們獲得的資訊會越多，這就是資訊含量 Entropy 的基本想法，Entropy 定義是某個隨機事件平均每個結果發生時所能傳達出的資訊量，公式如下

S = - \sum_{i = 1}^{n} P (x_{i}) l o g_{b} P (x_{i}) = \sum_{i = 1}^{n} P (x_{i}) I (x_{i})

Entropy 實際上是 Self-information

I (x_{i})

的期望值，Self-information 即每個隨機事件發生所傳達的資訊量，定義上資訊量與隨機事件發生的機率成反比，且機率爲 1 時資訊量爲 0，根據定義我們發現 log 有這樣的性質，因此

I (x_{i}) = l o g_{b} (\frac{1}{P (x_{i})})

帶入 P 值可以發現 P 越小 I 會越大，符合我們對資訊量的定義。

此外使用 log 還有一個性質：資訊是可加的，兩個獨立事件發生的資訊量為個別資訊量相加。

\begin{aligned} I (A + B) & = l o g_{b} (\frac{1}{P (A) P (B)}) \\ = l o g_{b} (\frac{1}{P (A)}) + l o g_{b} (\frac{1}{P (B)}) \\ = I (A) + I (B) \end{aligned}

觀察 Entropy 公式可以發現到，對於一分佈而言，因

\sum_{i = 1}^{n} P (x_{i}) = 1

，Entropy 的極大值發生在

P (x_{1}) = P (x_{2}) = . . . = P (x_{n}) = \frac{1}{n}

時

(S = l o g_{b} n)

，此分佈即爲 uniform distribution，因爲每個事件的發生的機率都一樣，不論發生哪個事件我們都可以獲得最大的資訊量，這就是資訊觀點的亂數，從剛剛文法的例子去想也很直覺，若我們難以從前面的資訊預測接下來會出現的詞，我們更有可能認為該句子是亂數產生的。

上面我們從兩個不同的觀點切入並獲得相同的結果，並有以下結論，隨機具有以下特性

每個隨機事件的發生機率相同
每個事件發生的機率互相獨立

並且我們可以用 Entropy 來量測亂度

S = - \sum_{i = 1}^{n} P (x_{i}) l o g_{b} P (x_{i})

PRNG

在現實世界中是否存在真正理想的亂數，目前仍有爭議，因為古典的科學建立在任何事件發生皆有因果關係之上，如果我們能夠獲取夠全面的資訊的話，就應該可以預測接下來會發生的事情，若存在一事件發生與過去所有狀態無關，則會打破這樣的因果關係；但近代的量子物理開始有發現到一些自然界的隨機性。

儘管如此我們可以透過一些手段，讓亂數雖然是經過一連串確定性的過程產生，卻可在我們使用的範圍內「不可預測」，這就是「偽亂數」(pseudo-random number)，產生偽亂數的產生器即為 pseudo-random number generator (PRNG)。

PRNG 的數學描述如下

考慮一函數集合

A = {A : {0, 1}^{n} \to {0, 1}^{*}}

為一系列檢測一組長度為 n 的序列是否為亂數的統計試驗，稱為 Adversaries，則有一函數

G : {0, 1}^{ℓ} \to {0, 1}^{n} with ℓ \leq n

稱為 Generator，試圖欺騙所有

A

裡面所有的

A

，並容許偏差

ϵ

，即對於所有

A

裡面的

A

，

A (G (U_{ℓ}))

和

A (U_{n})

的統計距離不大於

ϵ

，其中

U_{k}

為

{0, 1}^{k}

的 uniform distribution。

Adversaries 的例子可參考 Statistical randomness - Wikipedia#Tests

Randomness test

linux 底下有個命令 ent 可計算一輸入字串的 Entropy。

嘗試計算 linux 內建的 PRNG /dev/random

$ 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).

1 byte char 的大小為

2^{8}

，則最大的 entropy 為

S = l g (2^{8}) = 8

，可觀察到這次取出的結果相當接近這個理論值。

輸出的第二的部分為最佳壓縮率計算，這個值是從 entropy 計算而來的。

輸出的第三部分為卡方檢定的結果，卡方中的適合度檢定(Goodness of Fit)可用來檢定資料是否符合某種分佈，我們可以用卡方來檢驗我們的資料是否符合 Uniform distrubution，卡方檢定的 null hypothesis 如下

H_{0} : p_{i} = p_{i 0}, i = 1, . . ., k

H_{1} : 不 是 所 有 p_{i} 都 等 於 p_{i 0}

經由卡方檢定計算出的 p 值若小於設定的信心水準 (如0.05) 表示拒絕此虛無假說，則可宣稱此組資料不符合指定分佈，反之則宣稱無證據顯示此資料不符合該分佈。

當我們拒絕 null hypothesis 時我們明確表示
$H_{1}$ 爲真，但若檢定結果爲不顯著，結論並非「
$H_{0}$ 爲真」，而是「無證據顯示
$H_{1}$ 爲真」DanielChen

以這次執行的結果而言 p 值為 47.82%，離 5% 信心水準相當遙遠，因此判斷這次產生的亂數有可能符合 Uniform distribution。

最後三個部分分別為平均值檢定、蒙地卡羅法計算Pi值、及序列相關係數(可檢驗序列上每個值前後的相關性)，也是檢驗該序列是否是亂數的依據。

除了 ent 之外，還有一個很爲人所知的工具叫做 Diehard tests，他是一系列 Randomness test 的集合，最早由 George Marsaglia 開發，之後釋出的 linux 命令 dieharder 還有後續加入其他的亂數測試。

節錄部分執行結果

#=============================================================================#
#            dieharder version 3.31.1 Copyright 2003 Robert G. Brown          #
#=============================================================================#
   rng_name    |rands/second|   Seed   |
   /dev/urandom|  2.16e+07  | 670031799|
#=============================================================================#
        test_name   |ntup| tsamples |psamples|  p-value |Assessment
#=============================================================================#
   diehard_birthdays|   0|       100|     100|0.80323410|  PASSED  
      diehard_operm5|   0|   1000000|     100|0.99712783|   WEAK   
  diehard_rank_32x32|   0|     40000|     100|0.13150830|  PASSED  
    diehard_rank_6x8|   0|    100000|     100|0.60447291|  PASSED  
   diehard_bitstream|   0|   2097152|     100|0.14568798|  PASSED  
        diehard_opso|   0|   2097152|     100|0.49901549|  PASSED  
        diehard_oqso|   0|   2097152|     100|0.08646324|  PASSED  
         diehard_dna|   0|   2097152|     100|0.75489321|  PASSED  
diehard_count_1s_str|   0|    256000|     100|0.95499444|  PASSED  
diehard_count_1s_byt|   0|    256000|     100|0.59438482|  PASSED  
 diehard_parking_lot|   0|     12000|     100|0.00000000|  FAILED  
    diehard_2dsphere|   2|      8000|     100|0.00000000|  FAILED  
    diehard_3dsphere|   3|      4000|     100|0.00000000|  FAILED  
     diehard_squeeze|   0|    100000|     100|0.00000000|  FAILED  
        diehard_sums|   0|       100|     100|0.00000000|  FAILED  
        diehard_runs|   0|    100000|     100|0.96701019|  PASSED  
        diehard_runs|   0|    100000|     100|0.80232209|  PASSED  
...

Stat Trek RNG

以下用 JavaScript 簡單寫了一支小程式來計算 Random Number Generator 這個網站的 entropy 以及卡方值，新增一個書籤並貼上貼上以下代碼，名稱任意，進入網站並產生亂數表之後執行剛剛新增的書籤，結果會顯示在 browser console 裡面

javascript:(function(){var script=document.createElement('SCRIPT');script.src='https://rawgit.com/team6612/RandomnessTest/master/test-compress.js';document.body.appendChild(script);})()

執行結果(範圍0~99999,產生1000個)

Calculated entropy = 9.955784284662016
Ideal entropy for N=100000 is 16.609640474436812

Chi Square = 99999.99999987465

就 entropy 而言這次產生出來的亂數表距離理想值還有一段距離。

因為卡方分佈的 CDF 不太好實作，這邊只算了卡方值，可用查表的方式查出自由度 99999，

X^{2} = 99999.99999987465

時，

p = 0.49851 = 49.8 %

，

H_{0}

不顯著，無證據說此分佈不均勻。

但這邊有一個問題，取樣的數目可能太少而範圍太大，發生 local randomness 不足的可能很高，調整亂數範圍為 0~99

Calculated entropy = 6.559046277841441
Ideal entropy for N=100 is 6.643856189774724

Chi Square = 117.6

發現 entropy 與理想值相當接近，表示該亂數產生器產生出來的亂數有一定程度的亂度。卡方檢定部分，

p = 0.09787 = 9 %

，

H_{0}

不顯著，無證據說此分佈不均勻。因此結論此 RNG 通過 entropy 與卡方檢定，我們能在統計上確保一定程度的亂度。

Image Not Showing Possible Reasons
The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported
Learn More →

論文〈Dude, is my code constant time?〉重點提示

此論文提出一套量測方法與對應原始碼 dudect，用以量測程式執行時間是否能在常數時間

O (1)

完成，其特點如下

不侷限某特定 CPU 硬體平台，此特點相當重要。該論文也提及其它相關研究，往往只能在某特定 CPU 架構才能量測 (因為需要更新 CPU 的 microcode)。
該論文使用 Student's t-distribution (使用情境: 取樣大小不大且標準差無法得知的情況下) 之 Student's t-test 方法

利用 leakage detection test 測試程式是否是 constant time ，優點是在這個方法下不需要硬體模型，分為以下三個步驟進行

Repeatedly measure exeution time
分別用兩類(class) input data ，論文中提到的是 fix-vs-random ， fixed class 是為了 corner-case
Apply post-processing
- Cropping: 去掉某些極端值
- High-order preprocessing
Apply statistical test
- null hypothesis: "the two timing distributions are equal"
- if statistical test reject the null hypothesis, then we know it's not run in constant time
- Welch's t-test: 本專案使用的測試

Welch’s t-test

先簡介 Student’s t-distribution 就是抽樣下的常態分佈，抽樣越多就會越接近標準常態分佈(因為越接近母體)，圖中的 one-tailed test 和 two-tailed test 和 alternative hypothesis 有關，我們只需要看兩組數據是否不一樣，所以用 two-tailed test

alternative hypothesis: "the two timing distributions are not equal"

接著是 Welch’s t-test

先計算出兩個 class (class 0, class 1) 的平均值和變異數
Welch's t-test defines the statistic t by the following formula:
- $t = \frac{\bar{X_{0}} - \bar{X_{1}}}{\sqrt{\frac{V a r_{0}}{N_{0}} + \frac{V a r_{1}}{N_{1}}}}$
- $t$ 越接近 0 代表兩組數據的差異越小
- lab0 實作上只有到這裡，然後直接檢查
  $t$ 是否大於 t_threshold_moderate
approximate df (degree of freedom) using the Welch–Satterthwaite equation
- 由 df 我們可以得到相應的 t-distribution 如上圖的鐘型曲線
- 圖上的藍白交界處分別為 ±t
- 把藍色部份積分起來可以得到 p value 若 p 小於特定值 (通常 0.05 這部份沒有深入理解) 我們可以說兩組數據有統計上的顯著差異 (statistically significant difference)

實作細節

lab0-c 的 Welch's test 測試項目有兩個:

is_insert_tail_const(): 用來測試將字串加到 Queue 尾端的執行時間是否能在常數時間完成。
is_size_const(): 用來測試取得 Queue 大小時間是否能在常數時間完成。

每個測項又細分兩個階段:

Pre-processing: 主要測量待測物 (給定的程式碼) 的執行時間 (CPU cycle)，其步驟如下:
- 兩種產生 Queue 方式 (即對照 Welch's Test 所需兩組取樣資料):
  - Queue 大小 = 0
    $\to$ 論文提及的 fix class
  - Queue 大小 = n (隨機產生)，並插入 n 個相同的隨機字串至 Queue 開頭 (字串長度=7bytes, 加上結束字元 '\0' 則為 8bytes)
    $\to$ 論文提及的random class
- 執行待測物 (is_insert_tail_const() 或 is_size_const())，並記錄執行時間 (CPU cycles)
Statistical test: 將所記錄的CPU執行時間 (CPU cycles) 套用如下 Welch's test 公式，用以計算t value。可對照 lab0-c 原始碼的 t_push() 與 t_compute()。其所計算得出 t value 就可用來判斷測試項目是否在常數時間完成。

$t = \frac{\bar{X_{0}} - \bar{X_{1}}}{\sqrt{\frac{s_{1}^{2}}{N_{1}} + \frac{s_{2}^{2}}{N_{2}}}}$

$t$ 越接近 0 代表兩組數據的差異越小