執行人: kuanyu0712, Hlunlun
期末專題影片連結
rota1001
在 ldexpf 的實現中,需要去限定 exp
的範圍:
是否想過沒有分支的解決方案?像是用 -!!(new_exp >> 8)
可以生成出一個全為 1 或全為 0 的 bitmask,可以判斷是否 overflow,也可以用相似的方法去判斷 new_exp
的正負號。
另外是在浮點數運算部份要如何驗證計算結果是否與 IEEE 754 中定義的浮點數運算相同?(指的不是絕對值大小的關係,而是是否符合 IEEE 754 下的運算結果)。譬如說使用 math.h
中的函式去進行一些結果的比對,這是我在這個專題沒有看到的。然後在幾何平均數的部份,在專題解說裡面是說 可以忽略。從數值上來看是的,但是他是放在指數位,所以其實對於 mantissa 的影響蠻大的,這樣的計算方式會讓 mantissa 一定是 0。另外 e_sum /= n
這個部份做整數除法也會造成誤差。
EricccTaiwan
PELT 是一種策略,可以讓 scheduler 知道在 CPU 中執行的任務的重要性,這個重要性會用分數呈現…
關於 PELT,這裡所謂的「重要性」,具體是指 sched_entity
(se
) 中的哪個欄位或屬性?
eleanorLYJ
第二周的測驗 3 的範例寫得模糊。
以上例子應註明 custom_signal_handler 是單參數函式,若使用三個參數的自定義 handler function,則要寫入 sa_sigaction
,而非 sa.handler
裡,並且 sa_flags
中設定 SA_SIGINFO
。
weiso131
接下來關於兩個 32 bits 數值合在一起的部分,既然 y 在一開始就限定是 uint32_t 的型態,所以關於原先以下部分(-1L >> 32)去做 y 的遮罩我著實想不到什麼情形會發生從 uint32_t 轉型成 uint64_t 型態後,最高32個位元會發生非0的狀況
根據 C99 規範(ISO/IEC 9899:1999) §6.5.7「Shift operators」 第五項
The result of E1 >> E2 is E1 right-shifted E2 bit positions. If E1 has an unsigned type or if E1 has a signed type and a nonnegative value, the value of the result is the integral part of the quotient of E1 / 2E2 . If E1 has a signed type and a negative value, the resulting value is implementation-defined.
負值的右移結果可由不同編譯器/硬體平台自行決定,例如我在 x86 架構電腦上使用 gcc 13.3.0 ,得到的結果就是 0xffffffffffffffff
tsaiiuo
在 mpi_fdiv_qr
的範例中
在表示二進位的時候應該要新增 2 ,不然會讓範例程式碼閱讀起來會有些困惑,像是 11 ≥ 11 實際上應該改為 112 ≥ 112 或 310 ≥ 310,
MikazukiHikari
TODO: 第 2 週測驗題
應該改成第 3 週測驗題,畢竟連結和內容全都是第 3 週的。
wurrrrrrrrrr
關於第二週提到的 ceil_div
函數,沒有說明到為什麼在計算時,n 需要加上 d - 1 這個值,而不是直接加 d 或 d - 2 等其他值?
Mike1117
第一步是先做 alignment,原先我並不了解為什麼要限制低位元的 bits 為 0 才算是對齊,所以我參考了 你所不知道的 C 語言:記憶體管理、對齊及硬體特性 發現只要跨區取位元就會導致存取速度降低
未說明如何做 Alignment 。
Ian-Yen
在
mpi_compact
的程式碼中看起來沒有對rop->capacity == 0
的情況做處理,是否會導致導致 data[i] 越界?
HeLunWu0317
TODO: 第 3 週測驗題 測驗一
請問當 realloc 失敗,會被設為 NULL 後會做 absort ,那麼這裡不需要釋放記憶體嗎?
測驗題 1, 2 和延伸問題
不要列出參考題解,專注在程式碼本體。
以下程式碼是要可以用 32 bits 儲存的整數來達成任意長度整數的運算,所以會用 mpi_t
的結構來拼接成更多位元。
mpi_t
structure 示意圖:mpi_t
記憶體空間mpi_t
裡 data
所指向的所有記憶體mpi_t
裡位元空間不夠放時就會透過 capacity
指定的數量,擴展 data
的空間mpi_t
裡不需要的 data
空間釋放32 bits
的整數拆成多個 32 bits 單位存到 rop->data[n]
裡,給mpi使用關於記憶體重新配置可以看到以下 C99 的說明
C99 7.22.3.5 The realloc function
void *realloc(void *ptr, size_t size);
The realloc function deallocates the old object pointed to by ptr and returns a pointer to a new object that has the size specified by size.
因為 realloc 是以 bytes 為單位,所以下列這行就是因為一個 capacity
代表一個 uint32_t
的空間,裡面是存放 32 bits 的大小(= 4 bytes),所以才會是 capacity * 4
,因此我覺得乘以 sizeof(uint32_t)
或許會更直觀
接著再根據 C99 同段落的描述
The realloc function returns a pointer to the new object (which may have the same value as a pointer to the old object), or a null pointer if the new object has not been allocated.
可以發現若記憶體配置失敗會回傳 NULL
,因此可透過以下方式檢查,capacity != 0
代表有要求配置記憶體,但 !rop->data
仍然失敗,所以呼叫 abort()
終止程式
接著新分配出來的空間並不會有值,所以要主動初始化避免後續做運算時出錯。然後原本的空間因為 realloc
函式並不會更改到,會保留同一塊,所以不做操作
透過學習 mpi_enlarge
函式,mpi_compact
函式中一開始有以下這段的技巧看不懂,所以我先去理解 i != (size_t) -1
的部分
了解到 size_t
能表示的最大值就是 (size_t) -1
。因此 i=0 最後一次透過 i--
就會再繞回去變成 size_t
可表示的最大值,也就是 (size_t) -1
所以達成從 capacity 的 最大值 - 1
遍歷到 0
的效果。
接著只要 data[i]
沒有存放值就會把 capacity-1
以利後續收回空間
最後 assert()
的部分有參閱C99
C99 7.2.1.1 The assert macro
When it is executed, if expression is false, the assert macro writes information about the particular call that failed. It then calls the abort function.
若當中的表達式不成立主要會做兩件事,印出錯誤和呼叫 abort()
函式,所以就能防呆,確保 capacity
不會變成 (size_t)-1
這麼大的數值後再去呼叫realloc()
配置一個超大空間。
以 ceil_div(64, 31)
呼叫為例,64 bits 要拆成幾個 31 bits 單位的計算方式會是 (64 + 31 - 1) / 31
,因為透過C語言整數除以整數的特性會無條件捨去小數,因此會等於 3
為什麼不行
n / d + 1
因為當相除是整數時 反而會讓答案+1
參照 Linux 核心原始程式碼的整數除法
接著這部分就是把 64 bits
的 op
先取出最低位的 31 個位元放入data[n]
,所以INTMAX
的值應當就是 0x7FFFFFFF
作為遮罩,然後再往右移 31 bits 繼續就可以把 op
分成個別存放的 data[n]
m
組 乘以第 n
組對應的新位置 k
就會是 m + n
,如下圖說明接下來探討 carry 的部分:
tmp->data[k] += (r & INTMAX) + c
第 k 個位置就會是相乘完後結果 r 的最低 31 個位元tmp->data[k] &= INTMAX
確保值會落在31個bits以內理解到 carry 有值會發生在 tmp[k] 已超過 2 bits 時,因此可以擴展到 31 bits 的情況
TODO: 更詳細說明 mpi 除法和gcd
此為除法和餘數運算,參數定義為n = q × d + r
例子:
比如:
SWAR 把多個字元放進同一個暫存器,然後可以用一道指令同時處理,這題的話是一次比對多個位元是否等於目標字元,加速 memchr 搜尋
一開始有先提到兩個數字要判斷是不是奇數的作法,以下這個方法就是可以直接比對 64 bits 的數值
接下來關於兩個 32 bits 數值合在一起的部分,既然 y 在一開始就限定是 uint32_t
的型態,所以關於原先以下部分(-1L >> 32)
去做 y 的遮罩我著實想不到什麼情形會發生從 uint32_t
轉型成 uint64_t
型態後,最高32個位元會發生非0的狀況
所以我認為可以改成以下寫法,強制把 y 轉型成 uint64_t
,不用再額外對最高 32 bits 做處理
接下來把 SWAR 的技巧運用到原先是單個字元比對的 memchr 函式
void *memchr(const void s[.n], int c, size_t n);
The memchr() function scans the initial n bytes of the memory area pointed to by s for the first instance of c. Both c and the bytes of the memory area pointed to by s are interpreted as unsigned char.
找出在記憶體區段 s 裡面前 n 個 byte,第1個出現 byte 值為 c 的位置。
第一步是先做 alignment,原先我並不了解為什麼要限制低位元的 bits 為 0 才算是對齊,所以我參考了 你所不知道的 C 語言:記憶體管理、對齊及硬體特性 發現只要跨區取位元就會導致存取速度降低
接下來是要建立遮罩可以一次對整個 word 去做檢查
d
先往左 8 個bits 再 or 原本自己就會變成重複兩次BBBB
完全跟上面的操作是一樣的接著是要直接用 XOR 來比對是否有出現 0 (代表欲查找的c
字元有出現),但要如何快速知道一個 word 中是否有 bit 為 0,所以使用以下巨集,就是呼叫 DETECT_CHAR(X, mask)
時會直接被展開成 DETECT_NULL((X) ^ (mask)
然後又因為 DETECT_NULL
有被定義成對應 32 bits 或是 64 bits 的特殊處理(說明如 quiz3 中提及),最終有值代表找到了目標字元,零代表沒找到
在理解 signal(7) 的傳遞過程前,有先實際寫一段程式 (補github連結),設定按 ctrl+C
會讓行程跳去客製化的Signal Handler。
因此在這過程中我先認識到 Linux Signal 是什麼,參考 Linux Signal 解說
所以一開始我不懂 kernal space 究竟是在何時知道 signal disposition 被改變了,是等 signal 被送到 user space 後才開始處理嗎?或是是從 user space 透過 sigaction 告知 kernal space?
結論是透過 sigaction()
傳遞更改後的 struct sigaction
告知 kernal space 因為 sigaction()
是一種 system call,所以釐清傳遞方向後,以及設定 sigaction()
的前置作業,想繼續探討實際上運行時的流程
kernal space 在每次轉移回 user space 時,都會檢查是否有正在等待的 signal,
Pending Signals
從userspace送到kernal space? 不對, 是相反 比如ctrl+C 就是從kernal space 送停止的signal到正在執行的行程(程式?program?
收到signal後會有一些固定的反應
客製化收到signal後要幹嘛
寫在 Signal Handler 裡
Execution of signal handlers
會先把原本做的事暫存(在哪?
跑去剛剛寫的Signal Handler
把剛剛暫存的狀態復原
所以接下來怎麼「存檔」透過 getcontext() 會打包( 包含Program Counter Stack Pointer 和 CPU register的值?) 後推入 userspace stack 所以寫不好會造成stack overflow???
復原的時候就是用 setcontext() 去讀檔,跳回 Signal Handler 之前繼續做事
跟原先 setjmp()/longjmp() 舊的差異?
這段說明作業系統是如何控制收到 Signal 時的行為,並解釋了當你自訂處理流程時,如何處理暫停與恢復,確保程式順序(?不會錯亂
測驗題 和延伸問題
可以用來檢驗資訊傳遞完整與否的方式,是一種多項式的表示,計算方式採 mod 2 也就是可以透過 XOR 運算,所以加法和減法視為相同。mod 2 加法和減法運算定義為如下,所以才可以使用 XOR 運算:
操作的方法為:
發送端指定一個數值作為除數,並對訊息進行除法運算來取得餘數,然後將餘數附加到訊息的末尾。接收端對該訊息執行相同的運算來驗證是否有錯誤,如果餘數不為零,則表示訊息有誤。 節錄至 2025q1 第 4 週測驗題 所描述
不同於傳統 CRC-32,CRC-32C 主要是應用在硬體加速,此次主要使用反轉多項式 0x82F63B78
來實作 CRC-32C
類別 | 多項式 | 反轉多項式 | 應用 |
---|---|---|---|
CRC-32 | 0x04C11DB7 | 0xEDB88320 | ZIP, GZIP, PNG |
CRC-32C | 0x1EDC6F41 | 0x82F63B78 | SCTP,iSCSI, SSE4.2, |
由於我們最後的目的就是要得到 CRC 值,而這個過程是
所以我們需要進行除法,而這個過程是使用 XOR 來模擬 TODO: 為什麼XOR可以模擬除法?
同餘
以下舉例資料和多項式如何使用 XOR 模擬除法:
POLY
和data
的階會影響輸出?
可以參考 Cyclic Redundancy Check (CRC) Step-by-Step Calculator 逐步了解多項式運算過程,跟一般除法方式不一樣的地方是相減改為 XOR 運算,如第一步原本是 1101 - 1011 = 0010
,但這邊是做 XOR 1101 ^ 1011 = 0110
,最後得到的 011
即為 CRC
驗證: 將最後的餘數 011
加到資料 11010011
後方變成驗證所需被除數 11010011011
,除以相同除數 1011
,若是餘數是 0 就代表資料在傳遞過程中沒有損壞
至於為什麼將餘數附加到訊息的末尾後再用相同的除數會整除,參考以下數學式子。
所以把原多項式補0(也就是往左位移)可以表達成:
然後傳送的資料可以表達成:
接著把 代換掉:
最後可以發現多出兩個餘數相加,但根據 XOR 運算的特性
因此結論是若資料傳遞正確,必會整除。
(下方圖待改)
對於硬體來說,從 最高有效位元 MSB 做除法的方式相對於從 最低有效位元 LSB 複雜,所以我們需要對 POLY
做反轉進而從 LSB 開始做除法
反轉 POLY
成 1101
,如果最低位是 1 就進行 XOR,如果是 0 就只右移,重複 8 次,最後得到餘數 100
POLY
除數需要翻轉,原先的資料不需要,等全部做完後,要再把最終的 CRC 做反轉(0變1, 1變0)0x82f63b78
(是原先 0x1EDC6F41
的 bit reverse)參考 16 bit CRC Calculation | Error Detection| Modbus RTU | Using Reflected (or Reverse) Polynomial
傳統 CRC32 每次一個 bit 檢查太慢了,如果使用預先計算的 look up table ,可以一次處理多個位,從而進一步將速度提高
如果是 一次看 8 個位元 的資料做 XOR 的多項式 T
就會有 256 種組合,因為資料所有可能的 byte 值共有 種組合
所需的查找表總共要建立 256 個值
首先看到原本計算方式是判斷目前這個位元是否為 1,如果是代表可以用 XOR 消除(做除法),不是就直接向右移
等同於,也就是每次不是和 0 就是和 Polynomial
做 XOR
如果我們想要簡化這段程式碼可以有效運用 bitwise 運算,改成以下:
因為 crc & 1
只有 0 或 1 兩種可能,而 int
有 32 個位元,以下方式就可以有 32 位元的 mask 去取得當下是要和 Polynomial
還是 0
運算 XOR
而內部迴圈 for (unsigned int j = 0; j < 8; j++)
則是要模擬一次看 8 個位元的計算,雖然看起來運算量很龐大,但是這樣的建表只會進行一次,所以之後只要用查表就可以不用再計算一次,這樣就加快了之後 CRC 的運算
在查表時,可以一次除 8 個位元,再用這個值作為索引 index
去查表,查表在找的就是這個索引 index
迭代 8 次所得到的 crc 是多少,以下程式碼參考 Multiple Bits At Once
CRC VS. Fast CRC 過程以 POLY = 1101
舉例如下
雖然可以透過查表加速運算過程,但是這 256 個 32 位元的數值需要龐大的記憶體空間,所以這邊又進化到使用 Half-Byt Algorithm 節省使用記憶體空間
一次要看 8 個位元太耗記憶體,一次只看 4 個位元的查表就只有 16 個值,但是速度太慢,所以 Half-Byte 演算法採用了每個 byte 拆成兩個 4 個位元處理的設計,每次對一個 byte 的低 4 個位元和高 4 個位元各查表一次,這樣就和一次看 8 位元是一樣的了,對照以下示意圖可以看到 Half-byte Algorithm 兩次查表就對應傳統 CRC 在處理到第 4 位元和第 8 個位元時
data
就是一直往右移,所以可以做到一半再次查表因為一次只看 4 個位元,查表總共有 個值,和原本 256 個差了 16 倍,和原本建立方式一樣,只是每個值是處理 4 個位元的結果
所以想要用 Half-byte Algorithm 來實作 CRC-32C 的建表方式只要修改多項式 POLY
為 0x82f63b78
,並進行兩次查表即可
挑選 Linux 核心使用 CRC 的案例來說明 (執行 grep -r " CRC " *| less 命令),可說明以下 (用原始程式碼搭配)
輸入grep -r " CRC " *| less
命令後可以看到在 Linux 中尤其在 driver 有大量使用 CRC 的地方
(圖片來源)
CRC 就是落在 Data Link Layer確保實體傳輸的完整性
是 Intel 有線網卡的驅動程式(driver) ,就是在處理透過 PCIe 接上來的 Intel 網卡。當中的netdev.c
程式碼就是在處理:
Ethernet 協定有規定 CRC 的除數就是用以下多項式:
因此傳遞的過程是
CRCERRS
後續由 driver 來讀取來統計錯誤首先e1000_configure_rx
會根據封包大小來決定後續如何處理,有三種分別有三種
最常見的路徑會走e1000_clean_rx_irq()
,會先檢查 CRC 最後四個檢查碼是否有被刪除
開平方根 (square root, sqrt) 的整數運算、原理、加速機制,以及挑選 Linux 核心使用 CRC 的案例來說明,例如:
- drivers/net/wireless/ralink/rt2x00/rt2800lib.c
- drivers/video/fbdev/core/fbmon.c:
- drivers/md/bcache/bset.c:
透過 Linux 核心實例理解 CRC-CCITT(以 rt2800 驅動為例)
rt2800lib.c
是無線網卡的驅動程式,支援 Ralink RT2800 系列晶片
當中的 rt2800_check_firmware_crc(const u8 *data, const size_t len)
函式有用到crc_ccitt
驗證資料
首先在下列程式碼可以發現,最後兩個 bytes 為 CRC 驗證碼,所以先取出放入fw_crc
,因為該資料傳遞方式會 Big-Endian
參考 如何判斷Big-Endian or Little-Endian
然後參考 linux/lib/crc-ccitt.c 後,可以發現實際執行 crc 的是 crc_ccitt_byte
crc_ccitt(u16 crc, u8 const *buffer, size_t len)
再接著參考 linux-xlnx/include/linux/crc-ccitt.h中的crc_ccitt_byte(u16 crc, const u8 c)
可以發現在crc_ccitt_table[(crc ^ c) & 0xff]
是利用查表的方式計算 CRC,在 CRC-CCITT 裡,使用的 是 0x1021
,為 16-bit,但一次處理是 8 個位元,所以查表的索引會是 8 個位元,然後回傳的 CRC 會是 16 個位元,如下放上一小段crc_ccitt_table
的內容
所以最後再把crc = swab16(crc)
進行反轉再和原本的fw_crc
做比對
最後使用 swab16(crc) 把 16-bit 的 CRC 驗證碼進行 Byte 交換也就是 Endian 轉換,因為 Firmware 裡儲存的 CRC 與計算出來的 CRC Endian 順序不同
另外在以下的核心程式碼中有找到使用CRC的部份
drivers/usb/serial/safe_serial.c
在每個 USB 傳輸封包的末尾加上一段包含資料長度與 CRC10 驗證碼
一開始 程式碼 先定義 CRC
c
:新的 byte(fcs << 8) & 0x3ff
:左移 8 bit 並保留 10 bitscrc10_table[((fcs >> 2) & 0xff)]
:查表,index 是 fcs 的中間 8 bits^ (c)
:XOR 上新的輸入 byte在 static int safe_prepare_write_buffer(struct usb_serial_port *port, void *dest, size_t size)
中,是 CRC 附加的過程
要放 10-bit CRC,所以要把它:
步驟 | 說明 |
---|---|
count << 2 |
預留最後兩個 bytes,先放 count(協議需求) |
fcs_compute10() |
用 CRC10 算整包的檢查碼 |
fcs >> 8 |
高 2 bit,放進倒數第 2 byte 的低位 |
fcs & 0xff |
低 8 bit,放進倒數第 1 byte |
CRC 原理、如何加速 CRC 運算
測驗題 1, 2 和延伸問題
將程式碼大致分為三個區塊
acc = 0
: 所需遞迴深度已經完成,直接返回 0num >= base
: 代表 num
還可以繼續除以 base
,可以繼續遞迴num < base
: 計算num
的冪(指數為 base
),並繼續進行遞迴直到遞迴到達指定深度回傳對數,此程式碼對於輸入的限制是 base
和 num
都是大於 1 的任一浮點數
首先可以將計算 的部分使用快速冪計算,以下用 舉例,可以將指數部分 轉為二進位,再用指數的二進位展開來計算冪
實作如下,我們將指數 exp
當作二進位來看,跌代每個位元,如果該位元是 1,就可以將 num
對應的冪乘到 result
中,直到 exp
為 0
如此可以將原本計算 時間複雜度 減少到 ,整體計算對數的時間則是從 下降至
接下來為了讓任意正浮點數皆可以計算,這時候我們又想到了泰勒展開式,並且參考 fdlibm/e_log.c 講
關於 計算,將 拆解成以下
指數函數的計算 (特別是定點數)、原理,挑選 Linux 核心使用指數函數的案例,如 PELT 用到的 EWMA
PELT 是一種策略,可以讓 scheduler 知道在 CPU 中執行的任務的重要性,這個重要性會用分數呈現,而且這個分數會隨著時間用 指數加權移動平均 EWMA 衰減,也就是每個任務的重要性就會隨著時間下降。較常聽到的解釋是,PELT 是用來追蹤任務對 CPU 的負載計算策略,
Linux Kernel Document: Schedutil
With PELT we track some metrics across the various scheduler entities, from individual tasks to task-group slices to CPU runqueues. As the basis for this we use an Exponentially Weighted Moving Average (EWMA), each period (1024us) is decayed such that y^32 = 0.5. That is, the most recent 32ms contribute half, while the rest of history contribute the other half.
由於是每個週期是 毫秒,並且需要一個衰減因子 經過 次週期後剩一半
這個部分就是使用指數衰減的展現,以下說明指數加權移動平均 (EWMA)
Moving average 移動平均線 是指在 一整筆資料中選擇特定時間段的資料計算平均值
In statistics, a moving average is a calculation to analyze data points by creating a series of averages of different selections of the full data set.
EWMA 相較於一般應用於金融領域的 Simple moving average(SMA) 簡單移動平均線 不同的地方是在計算平均值時 SMA 每個資料權重是一樣的,都是 1
而 EWMA 會讓資料隨著時間權重降低,其中 smoothing factor 平滑因子 決定這個資料每次要衰減多少,且
exponential functions are used to assign exponentially decreasing weights over time
EWMA 會讓最新的資料 有較高的權重 ,其餘過去資料會因為不斷乘上 而變得越來越小,如下
隨著時間 可以忽略 這一項
這邊的 即可視為 PELT 所用到的衰減因子 , 的部分用 代替方便對應程式碼說明
在 Linux Kernel Document: Schedutil 文檔所提及的 EWMA 的計算將與幾何級數成正比的部分獨立出來成為一個函數
完整的 EWMA 就可以是 和 幾何級數 的比值,這樣也可以減少原本因為無窮級數的計算丟失精度
將分母與分子展開,可以得到我們一開始推導 EWMA 的公式,因為 ,所以分母 部分即為幾何級數,最後一定會收斂所以可以寫成
PELT 原始碼參考: linux/kernel/sched/pelt.c、pelt.h、sched-pelt.h、sched.h
相關解釋參考: Lin Yiwei 的 Linux 核心設計: Scheduler(4): PELT、蜗窝科技-PELT算法浅析
sched_avg
PELT 中主要用 EWMA 追蹤兩個指標: 任務的執行時間 running
、任務的等待時間 runnable
,這兩個指標的變化會反映在 CPU 被爭搶時
The Linux Kernel:PELT (Per Entity Load Tracking)
Using this(指EWMA) we track 2 key metrics: ‘running’ and ‘runnable’. ‘Running’ reflects the time an entity spends on the CPU, while ‘runnable’ reflects the time an entity spends on the runqueue. When there is only a single task these two metrics are the same, but once there is contention for the CPU ‘running’ will decrease to reflect the fraction of time each task spends on the CPU while ‘runnable’ will increase to reflect the amount of contention.
在 sched.h
中 sched_avg
用於定義一個計算負載平均的物件,其中跟 running 有關的數值就是 util_sum
(運用CPU的時間),與 runnable 有關的數值是 runnable_load_sum
(等待的時間),而負載是當前時間任務對系統產生的壓力,也就是 load_sum
應是 running + runnable 的時間
TODO: 這裡要計算的各種指標代表意義
decay_load
目的為計算
val 單位是?
也就是每個值每次衰減的計算, 代表第幾次衰減,我們期望最終在 負載平均週期數 LOAD_AVG_PERIOD
次衰減後剩原本的 倍
因為 kernel 中沒有 FPU , 直接算 較複雜,所以將 整合成以 為底數的冪,就可以用右移的位元操作完成這部分的指數運算
為了避免 因無法整除而丟失精度,首先觀察 ,商 為整數
利用指數性質
這樣就可以將 右移 ,並用查表方式來取得 ,實作可以看到 pelt.c#L49
其中 runnable_avg_yN_inv
即為查詢所有可能 的查表(將浮點數轉成 32 bit 的二進位結果),再用 mul_u64_u32_shr
將 和 這兩個值相乘取得最後衰減後的數值
在 include/linux/math64.h 有支持 128 位元整數和不支持的 mul_u64_u32_shr
兩種計算方式,如果沒有支持 128 位元整數的編譯器,可以用定點數的方式計算,實作如下
用定點數來解釋這段乘法,首先 a
拆解成高 32 位元的 ah
和 低 32 位元的 al
a
和 mul
的相乘即可分成高位元和低位元各自乘上 mul
後的和
先計算
再加上
shift
在這裡作用是甚麼 ? 如果看 shift
在整個衰減計算的位置,相當於
也就是可以防止衰減太快,有縮放衰減速度的作用
accumulate_sum
將所需計算的資料分為三個階段,d1
是上一個週期中剩下的時間,d1
是有經過完整週期的時間段的數量,d3
是目前週期已經過的時間,參考 Yiwei Lin 的 Linux 核心設計: Scheduler(4): PELT 說明 accumulate_sum
所以計算這三個階段總和 ,
加總所有時間到 ,包括過去時間片段 sa->period_contrib
其中 是目前已經過完的週期數量,一個週期是 微秒,如果已經過了 微秒,則 為
Step 1. 計算過去累積的負載 目前的衰減後的數值,即 ,程式碼如下,用 decay_load
更新累積至今的總和 *_sum
Step 2. 先更新 delta
成無法組成一個周期的剩餘時間,如果當前沒有負載的話可以不用計算 ,如果有則交由 _accumulate_pelt_segments
計算
最後加上新的負載貢獻,根據 ,最新的貢獻值不用乘上衰減因子,直接加上去就好
__accumulate_pelt_segments
計算當前負載貢獻值,即 d1
、 d2
和 d3
三個部分 ,也就是
d1
部分,即 ,有多少貢獻值就直接乘上衰減因子
d2
部分,即 ,由於這個計算會需要迴圈,這裡進行拆解以方便計算,
第一部分 就是無窮幾何級數,因為 最後會收斂
會是一個常數,所以定義此數 LOAD_AVG_MAX
在 linux/kernel/sched/sched-pelt.h,估算 如下
所以 大約是
有一些誤差
LOAD_AVG_MAX
是定義 47742
第二部分 展開如下
乘上 的結果,即使用 decay_load
來計算
第三部分即
d2
的計算完整如下
所以 d2
實作如下
d3
部分,最新進的貢獻,不用乘上衰減因子
最終回傳當前負載
___update_load_sum
___update_load_avg
以上函數都是和 相關的計算,完整的 最後還需除以
首先呼叫 get_pelt_divider()
, 相當於取得
最後相除,即完成 EWMA 的計算
get_pelt_divider
依照原先推倒的公式來說,應該是直接回傳
實作 ldexpf
根據 Floating-point API - The Linux Kernel documentation 中有提及:
Kernel code is normally prohibited from using floating-point (FP) registers or instructions, including the C float and double data types. This rule reduces system call overhead, because the kernel does not need to save and restore the userspace floating-point register state.
所以在 Linux 核心中可能無法使用 FPU 做運算,需要將浮點數轉成二進位制做運算
The function returns the result of multiplying the floating-point number x by 2 raised to the power exp.
先將 float
轉成 unsigned
然後再去取 sign
、exponent
、mantissa
根據 IEEE-754 Floating Point Converter 表示法去嘗試當 x
是 -1.2345 印出 sign
、exponent
、mantissa
但除了 sign
其他都不對,如下
檢討
關於浮點數和整數的轉換可以看到以下 C99 的說明,如果浮點數轉換成整數時會捨棄小數點的部分
C99 [6.3.1.4] Real floating and integer
When a finite value of real floating type is converted to an integer type other than _Bool, the fractional part is discarded (i.e., the value is truncated toward zero).
所以,當我用 cast operator 直接將 x
轉為 unsigned
時,x
就會只剩整數位的 ,雖然可以取到正確的 sign
,但根本取不到正確的 exponent
和 mantissa
改進本筆記的書寫,讓排版接近指定教材的風格
了解浮點數和整數轉換的問題後,我們知道不行直接將浮點數用 cast operator 轉換成整數,而是用 union
對同一塊記憶體取不同型態的數值,這種 type punning 類別重疊 的說明如下
C99 [6.5.2.3] Structure and union members 底下的補充說明
If the member used to read the contents of a union object is not the same as the member last used to store a value in the object, the appropriate part of the object representation of the value is reinterpreted as an object representation in the new type as described in 6.2.6 (a process sometimes called "type punning"). This might be a trap representation.
此處提及的 trap representation 在 C99 [6.2.6.1] 有解釋到:
問題
出自 C99 [6.2.6.1]
Certain object representations need not represent a value of the object type. If the stored value of an object has such a representation and is read by an lvalue expression that does not have character type, the behavior is undefined. If such a representation is produced by a side effect that modifies all or any part of the object by an lvalue expression that does not have character type, the behavior is undefined.41) Such a representation is called a trap representation
讀到這裡有個問題:那剛剛直接用 unsigned
去 cast 浮點數後輸出的 是代表 還是未定義 NaN
?
TODO: 所以是怎麼取值的
所以將浮點數轉為二進位數的方式如下
取出 sign
, exp
, frac
;
轉換為二進位置後,就可以來看要怎麼進行相乘運算
更多浮點數的運算可以看到教材 CS:APP 第 2 章重點提示和練習,所以取得 sign
、exp
、frac
就可以進行單精度浮點數和 的乘法運算了,以下使用 為例
首先浮點數編碼成二進位的方式為
也就是我們只要找到 (sign) 、 、 就可以運算浮點數和 的乘法, 已經知道是 ,將 轉換成二進位制會是無限循環,由於是單精度浮點數,所以就取到 位即可
Significand:
Exponent
Result
現在回到題目,我們這個時候會取得需要相乘的 2 的 次方,而 2 的次方 2 的冪 之間的相乘就是將次方數直接相加,如下
而此時的 會變成
那其實剛剛已經取出 sign
、exp
、frac
,所以可以直接將 n
加到 exp
那最後的答案就會是這三個值的 logical or
完全不會動到 frac
和 sign
,只要處理 exp
的部份,其實只要理解浮點數在電腦中的編碼模式就能解這題了
最後回傳 val.f
就是答案
CS:APP 第二章 有提到浮點數編碼後會有三種狀態
在運算 new_exp
之前可以先檢查 x
是否為有效值
x
是 NaN
,其實就是在檢查 exp
的部份
x
是無限大或無限小
x
是 0
再來需要檢查 new_exp
是否為有效值
new_exp
overflownew_exp
是否超過 也就是二進位制的
new_exp
underflow
x
是 normal 但輸出為 denormal
使用 lab0-c 規範的程式碼風格 (類似 Linux 核心原始程式碼)。
當你列出超連結時,應當給予對應的標題
只需檢查 exp
是否發生 overflow 或是 underflow ,當檢查到一位時就將 設置為 0,因為 任何數 與 正負無窮大 做運算皆會得到 NaN
(也就是未定義的數學結果)
要求:實作 geomean
,計算幾何平均數,不能用 FPU
之前已探討單精度浮點數的格式,此處可略過。
首先根據 IEEE 754
單精度浮點數佔用 32 位元的空間,其編碼方式為:
8.5 = 2³ + 1/2¹
轉成二進制為 1000.1
調整成 1.0001 x 2³
01111111 (127)
表示指數為0,所以指數為3時應先+3,表示成10000010 (130)
0001
因此放到 mantissa 裡就是00010000000000000000000
(共23bits)uint32_t a
的三個部分取出方法
將浮點數轉成定點數時,涉及數值誤差,應予以討論。
TOOD: 數值誤差
在 維基百科 介紹中描述了定點數表示法中,整數和小數的部分會用相同底數 來表示(就是會是一樣是十進位或是二進位),如果此有禮數 有理數的小數的部分到第 位,此有理數會是 的整數倍
In the fixed-point representation, the fraction is often expressed in the same number base as the integer part, but using negative powers of the base b. …Thus, if n fraction digits are stored, the value will always be an integer multiple of .
乍看之下,其實跟浮點數很相似,都是表示有理數的一種方式,但是可以看到,定點數會有 scaling factors 可讓數值表示成整數在運算,如下:
假設我們現在要將 和 相乘,我們可以想成這兩個定點數相乘
縮放因子在此次運算中即為
此時就可以不用用到 FPU 的將兩個整數相乘
最後再去乘上縮放因子即可得到結果
回到 Linux 核心中的浮點數相乘,我們如果可以用定點數的概念,那麼就可以用以下步驟得到結果
浮點數相乘
接著兩個單精度浮點數相乘
M_a = (1 << 23) | frac_a;
再做相乘注意用詞的精準,避免說「最前面」,而是「最高位元」
除了以上大概念外,還會遇到像是 mantissa 相乘時會有 overflow 和正規化的問題的問題
IEEE‑754 要求 mantissa 必須是 1.xxxx
的形式,也就是在 1.0~2.0 之間。如果乘出來 ≥ 2.0,就要把尾數右移一位,並把 exponent 加 1。
參考你所不知道的 C 語言: 浮點數運算中有提到浮點數的加法,若超過1也是把其中一個 mantissa 往右移,然後 exp 加一。
所以乘法應當也會有兩種情況
乘積大小 | Bit47 | 代表值範圍 | 處理方式 |
---|---|---|---|
≥ 2.0 | 1 | [2.0, 4.0) | P >>= 24; exp++ |
[1.0,2.0) | 0 | [1.0, 2.0) | P >>= 23; |
這裡都該是你的觀點,可省去說「我認為」。
減少項目縮排,讓書寫更流暢。
最後我認為 若要擴展到多個浮點數相乘,應該要一邊做正規化,保持 mantissa 在 23 bits,但這樣會一直遇到精準度遺失的問題。
TODO: 如何最大保留精準度
原本上面的作法是無條件捨取,所以誤差會很顯著,接著參考了 CMU: computer system 中 page 30 有提到可以採用 round 的方式,也就是對 mantissa
做四捨五入。
那實際上四捨五入有很多種實作的作法,根據 IEEE 754 和 你所不知道的 C 語言: 浮點數運算,有以下這幾種 mode,然後我們先選了 to nearest, ties from zero 去做四捨五入的實作。
Rounding mode | +11.5 | +12.5 | -11.5 | -12.5 |
---|---|---|---|---|
to nearest, ties to even | +12.0 | +12.0 | -12.0 | -12.0 |
to nearest, ties from zero | +12.0 | +13.0 | -12.0 | -13.0 |
toward 0 | +11.0 | +12.0 | -11.0 | -12.0 |
toward +inf | +12.0 | +13.0 | -11.0 | -12.0 |
toward -inf | +11.0 | +12.0 | -12.0 | -13.0 |
那 to nearest, ties from zero 可以想成理解成十進位的四捨五入也是看要 round 的後面那個位數而已。比如說一個 46 bits 要降成 23 bits,那就是在右邊數過來第 22 個 bits+1,因為這樣等同於對第 23 個 bit 做四捨五入,然後整體再往右移 23 個 bit ,程式碼如下: P = (P + (1ULL << (n - 1))) >> n;
注意書寫規範:中英文間用一個空白字元區隔。
以下圖片舉例 46 bits 欲捨棄 23 bits 的情況
參考: 留意浮點數運算的陷阱 和 to nearest, ties to even作法
將 全部相乘,並使用 取指數再進行開根號運算,如下:
初始化一個 ,假設 有 64 個位元
對傳入的陣列 進行迭代取出其數值的 部分,則其 的定點數為
用 會相乘每個 轉成定點數的數值部分
最後開根號
假設陣列
當
此時的 scaling factors
當
為了避免下一次和另一個 24 位元的數值相乘完溢位,需要將 47 位元進行捨入到 40 位元,這裡的捨入採用
此時的 scaling factors
當
此時的 scaling factors
最後開根號,因為浮點數的整數部分都是 1 ,所以 會是 0
似乎結果和第二個方法一樣,在沒有 FPU 的條件下可以只看 的部分就好
原先上面認為 要先全部相乘完再去做 取指數 ,但下面發現也能先對各別的 做 完再去相乘,因此我們想比較哪一個作法會保留更多的精準度。
如此,幾何平均數的運算可以變成 2 的冪的運算,所以在浮點數編碼成二進位後,就只要做指數部分的運算了
浮點數的 要轉換成底數為 2 的冪需要算出指數部分
Exponent 的總結果 會是
接著會發現 的部分計算結果恆為零,因為 本來就只有 23 個位元,乘上 就相當於右移 23 個位元直接變成 0,所以只剩 需要計算。
參考: What's the best way to calculate log(base2) of any number?
根據剛才的分析, log2
的部分就不用算了
第一個和第二個方法都沒有辦法到很精確的數值,因為 的部份最後都會被捨棄,即使取 也沒有辦法到很精確,所以可以用比較直觀的第二個方法實做
誤差分析
TODO: 如何比較兩種浮點數操作方法的精準度損失?
為了比較兩種實作的精準度,我們使用 double 作為 ground truth,接著比較兩種作法它們和 ground truth 的誤差,所以使用相對誤差作為指標:error = |approx - truth| / |truth|
expm1(x)
討論此處討論 linux kernel 中浮點數運算對 expm1()
的影響
exp(x)
要使用單精度浮點數還是雙精度浮點數進行分析?
在 IEEE 754 中對浮點數進行二進位編碼,所以可以將 轉為與 2 的冪相關的數就能更容易計算,參考 ulysses-libc/src/math/exp.c 如下
我們期望 是 倍的 並且加上誤差 ,所以將 x 分解成
所以 應近似於
則可以分解成
是 和 的誤差,也就是一個正確值 與近似值 的差
如此可以專注於需要高精度計算的 部分
解 和解 的數值範圍並不一樣,這種拆解方式可以防止直接算 時益位
以往計算 都是直接使用泰勒展開式
但是電腦沒有無窮大,所以會丟失精度,在參考的程式碼中定義了一個 special rational 函數
將 代入
由幾何級數的概念可推導出:
所以 前 5 項是
此時發現使用雷米茲演算法的五階多項式就能將誤差控制在 內
,,而其中的 5 個係數就是
P5 好像不是 1/30,240,000
由 得到
用 代入計算
此時討論如果 接近 0:
計算 時,分母中較大的 和較小的 相減可能會失去精度,如:
假設支持 15 位的精準度, :
但是因為只能支持 15 位,後面的位數直接被捨棄剩
所以整理成
因為 算起來會很小,所以 相對 會更接近 2,減少精度損失,而且可以讓 單獨計算
算出 後,就可以用 my_ldexpf
算出
expm1