Try   HackMD

Linux 核心專題: bitops 相關測驗題

執行人: kuanyu0712, Hlunlun
期末專題影片連結

Reviewed by rota1001

ldexpf 的實現中,需要去限定 exp 的範圍:

    if(new_exp >= 0xff){ // overflow or nan
        val.u32 = sign | 0x7f800000;
        return val.f;
    }
    
    if(new_exp <=0){ // underflow
        val.u32 = sign;
        return val.f;
    }

是否想過沒有分支的解決方案?像是用 -!!(new_exp >> 8) 可以生成出一個全為 1 或全為 0 的 bitmask,可以判斷是否 overflow,也可以用相似的方法去判斷 new_exp 的正負號。

另外是在浮點數運算部份要如何驗證計算結果是否與 IEEE 754 中定義的浮點數運算相同?(指的不是絕對值大小的關係,而是是否符合 IEEE 754 下的運算結果)。譬如說使用 math.h 中的函式去進行一些結果的比對,這是我在這個專題沒有看到的。然後在幾何平均數的部份,在專題解說裡面是說

log2(frac×223+1) 可以忽略。從數值上來看是的,但是他是放在指數位,所以其實對於 mantissa 的影響蠻大的,這樣的計算方式會讓 mantissa 一定是 0。另外 e_sum /= n 這個部份做整數除法也會造成誤差。

Reviewed by EricccTaiwan

PELT 是一種策略,可以讓 scheduler 知道在 CPU 中執行的任務的重要性,這個重要性會用分數呈現

關於 PELT,這裡所謂的「重要性」,具體是指 sched_entity (se) 中的哪個欄位或屬性?

Reviewed by eleanorLYJ

第二周的測驗 3 的範例寫得模糊。

struct sigaction sa;
sa.sa_handler = custom_signal_handler; // 指定要執行的函式
sa.sa_flags = 0;
sigaction(SIGINT, &sa, NULL); // system call

以上例子應註明 custom_signal_handler 是單參數函式,若使用三個參數的自定義 handler function,則要寫入 sa_sigaction,而非 sa.handler 裡,並且 sa_flags 中設定 SA_SIGINFO

Reviewed by 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

Reviewed by tsaiiuo

mpi_fdiv_qr 的範例中

0100  ← 商 (4)
   ┌────
111101  ← 被除數 (13)11↓↓  ← 1111,可以除,商寫1,餘數=000↓  ← 帶下來0,得00 < 11,商寫001  ← 帶下來1,得01 < 11,商寫01  ← 餘數 (1)

在表示二進位的時候應該要新增 2 ,不然會讓範例程式碼閱讀起來會有些困惑,像是 11 ≥ 11 實際上應該改為 112 ≥ 112 或 310 ≥ 310

Reviewed by MikazukiHikari

TODO: 第 2 週測驗題

應該改成第 3 週測驗題,畢竟連結和內容全都是第 3 週的。

Reviewed by wurrrrrrrrrr

關於第二週提到的 ceil_div 函數,沒有說明到為什麼在計算時,n 需要加上 d - 1 這個值,而不是直接加 d 或 d - 2 等其他值?

Reviewed by Mike1117

第一步是先做 alignment,原先我並不了解為什麼要限制低位元的 bits 為 0 才算是對齊,所以我參考了 你所不知道的 C 語言:記憶體管理、對齊及硬體特性 發現只要跨區取位元就會導致存取速度降低

未說明如何做 Alignment 。

Reviewed by Ian-Yen

mpi_compact 的程式碼中看起來沒有對 rop->capacity == 0 的情況做處理,是否會導致導致 data[i] 越界?

Reviewed by HeLunWu0317

TODO: 第 3 週測驗題 測驗一

請問當 realloc 失敗,會被設為 NULL 後會做 absort ,那麼這裡不需要釋放記憶體嗎?

TODO: 第 3 週測驗題

測驗題 1, 2 和延伸問題

測驗一整理

不要列出參考題解,專注在程式碼本體。

以下程式碼是要可以用 32 bits 儲存的整數來達成任意長度整數的運算,所以會用 mpi_t 的結構來拼接成更多位元。

/* mpi: Multi-Precision Integers */
typedef struct {
    uint32_t *data;                           
    size_t capacity;
} mpi_t[1];
  • mpi_t structure 示意圖:






G


cluster_mpi

mpi_t x



xdata

x.data → data array



d0

data[0]:32 bits



xdata->d0





d1

data[1]:32 bits



xdata->d1





d2

data[2]:32 bits



xdata->d2





xcap

x.capacity = 3



先整理各個函式的功用

  • mpi_init: 初始化 mpi_t 記憶體空間
  • mpi_clear: 在整個運算後釋放 mpi_tdata 所指向的所有記憶體
  • mpi_enlarge: mpi_t 裡位元空間不夠放時就會透過 capacity 指定的數量,擴展 data 的空間
  • mpi_compact: 把 mpi_t 裡不需要的 data 空間釋放
  • mpi_set_u64: 設定mpi,把超過 32 bits 的整數拆成多個 32 bits 單位存到 rop->data[n] 裡,給mpi使用

mpi_enlarge

關於記憶體重新配置可以看到以下 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) 或許會更直觀

- rop->data = realloc(rop->data, capacity * 4);
+ rop->data = realloc(rop->data, capacity * 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() 終止程式

if (!rop->data && capacity != 0) {
    fprintf(stderr, "Out of memory (%zu words requested)\n", capacity);
    abort();
}

接著新分配出來的空間並不會有值,所以要主動初始化避免後續做運算時出錯。然後原本的空間因為 realloc 函式並不會更改到,會保留同一塊,所以不做操作

for (size_t n = min; n < capacity; ++n)
    rop->data[n] = 0;

mpi_compact

透過學習 mpi_enlarge 函式,mpi_compact 函式中一開始有以下這段的技巧看不懂,所以我先去理解 i != (size_t) -1 的部分

(size_t)-1 
= UINT32_MAX 
= 2^32 - 1 
= 4294967295

了解到 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()配置一個超大空間。

for (size_t i = rop->capacity - 1; i != (size_t) -1; --i) {
    if (rop->data[i])
        break;
    capacity--;
}
assert(capacity != (size_t) -1);

ceil_div

ceil_div(64, 31) 呼叫為例,64 bits 要拆成幾個 31 bits 單位的計算方式會是 (64 + 31 - 1) / 31 ,因為透過C語言整數除以整數的特性會無條件捨去小數,因此會等於 3

static size_t ceil_div(size_t n, size_t d)
{
    return (n + d - 1) / d;
}

為什麼不行 n / d + 1
因為當相除是整數時 反而會讓答案+1
參照 Linux 核心原始程式碼的整數除法

mpi_set_u64

接著這部分就是把 64 bitsop 先取出最低位的 31 個位元放入data[n],所以INTMAX的值應當就是 0x7FFFFFFF 作為遮罩,然後再往右移 31 bits 繼續就可以把 op 分成個別存放的 data[n]

for (size_t n = 0; n < capacity; ++n) {
    rop->data[n] = op & INTMAX;
    op >>= 31;
}

mpi_mul_naive

for (size_t n = 0; n < op1->capacity; ++n) { 
    for (size_t m = 0; m < op2->capacity; ++m) {
        uint64_t r = (uint64_t) op1->data[n] * op2->data[m];
        uint64_t c = 0;
        for (size_t k = n + m; c || r; ++k) {
            if (k >= tmp->capacity)
                mpi_enlarge(tmp, tmp->capacity + 1);
            tmp->data[k] += (r & INTMAX) + c;
            r >>= 31;
            c = tmp->data[k] >> 31;
            tmp->data[k] &= INTMAX;
        }
    }                                     
}
  • 就像一般的直式乘法,每一位元都去乘對方的每一位元
  • 根據一般直式,確實第 m 組 乘以第 n 組對應的新位置 k 就會是 m + n,如下圖說明
           A2      A1      A0     (A 的位元)

×          B2      B1      B0     (B 的位元)
---------------------------------------------------
bit[3]  bit[2]  bit[1]  bit[0] 

  A0×B0  → bit[0]    (0+0)
  A0×B1  → bit[1]    (0+1)
  A0×B2  → bit[2]    (0+2)

  A1×B0  → bit[1]    (1+0)
  A1×B1  → bit[2]    (1+1)
  A1×B2  → bit[3]    (1+2)

接下來探討 carry 的部分:

  • 例子1 (該例子不會觸發c!=0的情況)
    • 當第一次時 tmp->data[k] += (r & INTMAX) + c 第 k 個位置就會是相乘完後結果 r 的最低 31 個位元
    • 假設前面還有 5 個位元,則r>>31後會只剩前面那五個
    • 接著c依然保持0
    • 然後tmp->data[k] &= INTMAX 確保值會落在31個bits以內
    • 接下來 r 有值會繼續第三層迴圈(r 或 c 還沒處理完就會繼續第三層迴圈)
  • 例子2
    就是把全部的數值用最大值,假設以下例子為每格只能裝 2 bits,對應到下列code
tmp[k] += (r & INTMAX) + c;
r >>= 2;
c = tmp[k] >> 2;
tmp[k] &= INTMAX;
op1 = [0b11, 0b11]
op2 = [0b11, 0b11]

🔸 op1[0] * op2[0] = 0b11 * 0b11 = 0b1001
  ➤ 位置 0:                        ➤ 位置 1:
    原本 tmp[0] = 0b00                原本 tmp[1] = 0b00
    r & INTMAX   = 0b01               r & INTMAX   = 0b10
    carry(c)     = 0b0000             carry(c)     = 0b0000
    sum          = 0b0001             sum          = 0b0010
    ➤ 新的 tmp[0] = 0b01              ➤ 新的 tmp[1] = 0b10
    ➤ 新的 carry   = 0b0000           ➤ 新的 carry   = 0b0000
    ➤ 剩餘 r       = 0b0010           ➤ 剩餘 r       = 0b0000

🔸 op1[0] * op2[1] = 0b11 * 0b11 = 0b1001
  ➤ 位置 1:                        ➤ 位置 2:
    原本 tmp[1] = 0b10                原本 tmp[2] = 0b00
    r & INTMAX   = 0b01               r & INTMAX   = 0b10
    carry(c)     = 0b0000             carry(c)     = 0b0000
    sum          = 0b0011             sum          = 0b0010
    ➤ 新的 tmp[1] = 0b11              ➤ 新的 tmp[2] = 0b10
    ➤ 新的 carry   = 0b0000           ➤ 新的 carry   = 0b0000
    ➤ 剩餘 r       = 0b0010           ➤ 剩餘 r       = 0b0000

🔸 op1[1] * op2[0] = 0b11 * 0b11 = 0b1001
  ➤ 位置 1:                        ➤ 位置 2:
    原本 tmp[1] = 0b11                原本 tmp[2] = 0b10
    r & INTMAX   = 0b01               r & INTMAX   = 0b10
    carry(c)     = 0b0000             carry(c)     = 0b0001
    sum          = 0b0100             sum          = 0b0101
    ➤ 新的 tmp[1] = 0b00              ➤ 新的 tmp[2] = 0b01
    ➤ 新的 carry   = 0b0001           ➤ 新的 carry   = 0b0001
    ➤ 剩餘 r       = 0b0010           ➤ 剩餘 r       = 0b0000

✅ 最終結果 tmp[]:
  tmp[0] = 0b01    tmp[1] = 0b00    tmp[2] = 0b10    tmp[3] = 0b11

理解到 carry 有值會發生在 tmp[k] 已超過 2 bits 時,因此可以擴展到 31 bits 的情況

TODO: 更詳細說明 mpi 除法和gcd

mpi_fdiv_qr

此為除法和餘數運算,參數定義為n = q × d + r
例子:

0100  ← 商 (4)
   ┌────
111101  ← 被除數 (13)
   │11↓↓  ← 1111,可以除,商寫1,餘數=0
   │ 00↓  ← 帶下來0,得00 < 11,商寫0
   │  01  ← 帶下來1,得01 < 11,商寫01  ← 餘數 (1)

mpi_gcd

mpi_fdiv_qr(q, r, op1, op2);  // op1 = q * op2 + r
mpi_gcd(rop, op2, r);         // 遞迴:gcd(op2, r)

比如:

gcd(48, 18) → gcd(18, 12) → gcd(12, 6) → gcd(6, 0) →  6

測驗 2 整理

補足背景知識

SWAR 把多個字元放進同一個暫存器,然後可以用一道指令同時處理,這題的話是一次比對多個位元是否等於目標字元,加速 memchr 搜尋

一開始有先提到兩個數字要判斷是不是奇數的作法,以下這個方法就是可以直接比對 64 bits 的數值

static uint64_t SWAR_ODD_MASK = (1L << 32) + 1;
bool both_odd_swar(uint64_t xy) {
    return (xy & SWAR_ODD_MASK) == SWAR_ODD_MASK;
}

接下來關於兩個 32 bits 數值合在一起的部分,既然 y 在一開始就限定是 uint32_t 的型態,所以關於原先以下部分(-1L >> 32)去做 y 的遮罩我著實想不到什麼情形會發生從 uint32_t 轉型成 uint64_t 型態後,最高32個位元會發生非0的狀況

 -1L → 0xFFFFFFFFFFFFFFFF
 (-1L >> 32) → 0x00000000FFFFFFFF

所以我認為可以改成以下寫法,強制把 y 轉型成 uint64_t,不用再額外對最高 32 bits 做處理

static inline uint64_t bit_compound(uint32_t x, uint32_t y)
{
-    return ((0L + x) << 32) | ((y + 0L) & (-1L >> 32));
+    return ((uint64_t)x << 32) | (uint64_t)y;

memchr

接下來把 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 語言:記憶體管理、對齊及硬體特性 發現只要跨區取位元就會導致存取速度降低

#define UNALIGNED(X) ((long) X & (sizeof(long) - 1))

接下來是要建立遮罩可以一次對整個 word 去做檢查

  • 首先先把要查找的字元 d 先往左 8 個bits 再 or 原本自己就會變成重複兩次
  • 然後接著再往左 16 bits 後 or 自己就會再重複兩次,也就是變四倍,這樣就足夠 32 bits 了
  • 所以接下來原先填空的 BBBB 完全跟上面的操作是一樣的
unsigned long mask = d << 8 | d;
mask |= mask << 16;
for (unsigned int i = 32; i < LBLOCKSIZE * 8; i <<= 1)
    mask |= mask << i; //BBBB

接著是要直接用 XOR 來比對是否有出現 0 (代表欲查找的c字元有出現),但要如何快速知道一個 word 中是否有 bit 為 0,所以使用以下巨集,就是呼叫 DETECT_CHAR(X, mask) 時會直接被展開成 DETECT_NULL((X) ^ (mask)

if (DETECT_CHAR(CCCC, mask))
#define DETECT_CHAR(X, mask) DETECT_NULL((X) ^ (mask))

然後又因為 DETECT_NULL 有被定義成對應 32 bits 或是 64 bits 的特殊處理(說明如 quiz3 中提及),最終有值代表找到了目標字元,零代表沒找到

#define DETECT_NULL(X) (((X) - (0x01010101)) 
& ~(X) & (0x80808080))
#define DETECT_NULL(X) \
    (((X) - (0x0101010101010101)) & ~(X) & (0x8080808080808080))

測驗 3 整理

補足背景知識

在理解 signal(7) 的傳遞過程前,有先實際寫一段程式 (補github連結),設定按 ctrl+C 會讓行程跳去客製化的Signal Handler。
因此在這過程中我先認識到 Linux Signal 是什麼,參考 Linux Signal 解說

  • Signal disposition: 代表 user space 中的 process 收到來自 kernal space 的 signal 會有幾個固定的行為反應(如終止行程、忽略該 signal或暫停該行程等等)
  • 如果想改變這個固定的反應可以透過 sigaction(2) 去設定,是一種system call ,有以下幾種行為可以選擇
    • 維持預設行為
    • 忽略該 signal
    • 透過 signal handler (programmer-defined function) 來處理該 signal
struct sigaction sa;
sa.sa_handler = custom_signal_handler; // 指定要執行的函式
sa.sa_flags = 0;
sigaction(SIGINT, &sa, NULL); // system call

所以一開始我不懂 kernal space 究竟是在何時知道 signal disposition 被改變了,是等 signal 被送到 user space 後才開始處理嗎?或是是從 user space 透過 sigaction 告知 kernal space?
結論是透過 sigaction() 傳遞更改後的 struct sigaction 告知 kernal space 因為 sigaction() 是一種 system call,所以釐清傳遞方向後,以及設定 sigaction() 的前置作業,想繼續探討實際上運行時的流程

Execution of signal handlers

kernal space 在每次轉移回 user space 時,都會檢查是否有正在等待的 signal,

Pending Signals

混亂的思緒整理中。

signal(7)

從userspace送到kernal space? 不對, 是相反 比如ctrl+C 就是從kernal space 送停止的signal到正在執行的行程(程式?program?

收到signal後會有一些固定的反應

sigaction(2)

客製化收到signal後要幹嘛
寫在 Signal Handler 裡

Execution of signal handlers
會先把原本做的事暫存(在哪?
跑去剛剛寫的Signal Handler
把剛剛暫存的狀態復原

userspace context

所以接下來怎麼「存檔」透過 getcontext() 會打包( 包含Program Counter Stack Pointer 和 CPU register的值?) 後推入 userspace stack 所以寫不好會造成stack overflow???
復原的時候就是用 setcontext() 去讀檔,跳回 Signal Handler 之前繼續做事
跟原先 setjmp()/longjmp() 舊的差異?

小結

這段說明作業系統是如何控制收到 Signal 時的行為,並解釋了當你自訂處理流程時,如何處理暫停與恢復,確保程式順序(?不會錯亂

TODO: 第 4 週測驗題

測驗題 和延伸問題

CRC (Cyclic Redundancy Check)

可以用來檢驗資訊傳遞完整與否的方式,是一種多項式的表示,計算方式採 mod 2 也就是可以透過 XOR 運算,所以加法和減法視為相同。mod 2 加法和減法運算定義為如下,所以才可以使用 XOR 運算:

1 + 1 = 0     0 - 0 = 0
1 + 0 = 1     0 - 1 = 1
0 + 1 = 1     1 - 0 = 1
0 + 0 = 0     1 - 1 = 0

操作的方法為:

發送端指定一個數值作為除數,並對訊息進行除法運算來取得餘數,然後將餘數附加到訊息的末尾。接收端對該訊息執行相同的運算來驗證是否有錯誤,如果餘數不為零,則表示訊息有誤。 節錄至 2025q1 第 4 週測驗題 所描述

CRC-32C(Castagnoli) 實作分析

CRC-32 VS. CRC-32C

不同於傳統 CRC-32,CRC-32C 主要是應用在硬體加速,此次主要使用反轉多項式 0x82F63B78 來實作 CRC-32C

類別 多項式 反轉多項式 應用
CRC-32 0x04C11DB7 0xEDB88320 ZIP, GZIP, PNG
CRC-32C 0x1EDC6F41 0x82F63B78 SCTP,iSCSI, SSE4.2,

XOR 模擬除法

由於我們最後的目的就是要得到 CRC 值,而這個過程是

data÷0x1EDC6F41=+CRC

所以我們需要進行除法,而這個過程是使用 XOR 來模擬 TODO: 為什麼XOR可以模擬除法?
同餘

以下舉例資料和多項式如何使用 XOR 模擬除法:

data: 11010011
POLY: 1011

POLYdata 的階會影響輸出?

可以參考 Cyclic Redundancy Check (CRC) Step-by-Step Calculator 逐步了解多項式運算過程,跟一般除法方式不一樣的地方是相減改為 XOR 運算,如第一步原本是 1101 - 1011 = 0010 ,但這邊是做 XOR 1101 ^ 1011 = 0110,最後得到的 011 即為 CRC

          11110001
     -------------
1011 ) 11010011000
       1011
       ----
        1100011000
        1011
        ----
         111011000
         1011
         ----
          10111000
          1011
          ----
              1000
              1011
              ----
               011           

驗證: 將最後的餘數 011 加到資料 11010011 後方變成驗證所需被除數 11010011011 ,除以相同除數 1011 ,若是餘數是 0 就代表資料在傳遞過程中沒有損壞


     -------------
1011 ) 11010011011
       1011
       ----
        1100
        1011
        ----
         1110
         1011
         ----
          1011
          1011
          ----
              1011
              1011
              ----
                 0

至於為什麼將餘數附加到訊息的末尾後再用相同的除數會整除,參考以下數學式子。

  • M(x)
    :原始資料的多項式
  • P(𝑥)
    :生成多項式(CRC 的除數)
  • Q(x)
    :商
  • xr
    :就是把 𝑀(𝑥)「後面補 0」

所以把原多項式補0(也就是往左位移)可以表達成:

M(x)xr=Q(x)P(x)+R(x)
然後傳送的資料可以表達成:
T(x)=M(x)xr+R(x)

接著把
M(x)xr
代換掉:
T(x)=Q(x)P(x)+R(x)+R(x)

最後可以發現多出兩個餘數相加,但根據 XOR 運算的特性
R(x)+R(x)=0

因此結論是若資料傳遞正確,必會整除
(下方圖待改)
image

對於硬體來說,從 最高有效位元 MSB 做除法的方式相對於從 最低有效位元 LSB 複雜,所以我們需要對 POLY 做反轉進而從 LSB 開始做除法

反轉位元組順序 Swap Endianness

uint32_t crc32_naive(uint32_t crc, uint8_t v)
{
    crc ^= v; // crc = crc ^ v;
    for (int bit = 0; bit < 8; bit++) {
        if (crc & 1)
            crc = (crc >> 1) ^ UINT32_C(0x82f63b78); // reversed_POLY
        else
            crc = (crc >> 1);
    }
    return crc;
}

反轉 POLY1101,如果最低位是 1 就進行 XOR,如果是 0 就只右移,重複 8 次,最後得到餘數 100

  • 只有POLY除數需要翻轉,原先的資料不需要,等全部做完後,要再把最終的 CRC 做反轉(0變1, 1變0)
  • 0x82f63b78(是原先 0x1EDC6F41 的 bit reverse)
  • 使用該方法的原因就是因為硬體是 LSB
    TODO:為什麼要先右移在 XOR、餘數不一樣?
                   crc          i
init CRC        00000000           
data           ^11010011           
                --------           
update CRC      11010011
shift           01101001
POLY           ^    1101
                --------
update CRC      01100100        0
shift           00110010        1
shift           00011001        2
                         
shift           00001100              
POLY           ^    1101               
                --------               
update CRC      00000001        3       
shift           00000000               
POLY           ^    1101               
                --------               
update CRC      00001101        4       
shift           00000110                 
POLY           ^    1101
                --------
update CRC      00001011        5
shift           00000101
POLY           ^    1101
                --------
update CRC      00001000        6
shift           00000100        7            

參考 16 bit CRC Calculation | Error Detection| Modbus RTU | Using Reflected (or Reverse) Polynomial

Fast CRC32

傳統 CRC32 每次一個 bit 檢查太慢了,如果使用預先計算的 look up table ,可以一次處理多個位,從而進一步將速度提高

如果是 一次看 8 個位元 的資料做 XOR 的多項式 T 就會有 256 種組合,因為資料所有可能的 byte 值共有

28=256 種組合

(crc >> 8) ^ T

所需的查找表總共要建立 256 個值

for (unsigned int i = 0; i <= 0xFF; i++)
{
    uint32_t crc = i;
    for (unsigned int j = 0; j < 8; j++)
        crc = (crc >> 1) ^ (-int(crc & 1) & Polynomial);
    Crc32Lookup[i] = crc;
}

首先看到原本計算方式是判斷目前這個位元是否為 1,如果是代表可以用 XOR 消除(做除法),不是就直接向右移

crc = (crc & 1) ? ((crc >> 1) ^ Polynomial) : (crc >> 1);

等同於,也就是每次不是和 0 就是和 Polynomial 做 XOR

crc = (crc & 1) ? ((crc >> 1) ^ Polynomial) : (crc >> 1) ^ 0;

如果我們想要簡化這段程式碼可以有效運用 bitwise 運算,改成以下:

crc = (crc >> 1) ^ (-int(crc & 1) & Polynomial);

因為 crc & 1 只有 0 或 1 兩種可能,而 int 有 32 個位元,以下方式就可以有 32 位元的 mask 去取得當下是要和 Polynomial 還是 0 運算 XOR

-int(1) = 0xFFFFFFFF;  ->  0xFFFFFFFF & Polynomial = Polynomial
-int(0) = 0x00000000;  ->  0x00000000 & Polynomial = 0 

而內部迴圈 for (unsigned int j = 0; j < 8; j++) 則是要模擬一次看 8 個位元的計算,雖然看起來運算量很龐大,但是這樣的建表只會進行一次,所以之後只要用查表就可以不用再計算一次,這樣就加快了之後 CRC 的運算

在查表時,可以一次除 8 個位元,再用這個值作為索引 index 去查表,查表在找的就是這個索引 index 迭代 8 次所得到的 crc 是多少,以下程式碼參考 Multiple Bits At Once

uint32_t crc32_1byte(uint8_t data,  uint32_t previousCrc32 = 0)
{
    uint32_t crc = ~previousCrc32;    
    uint32_t index = crc ^ data;
    crc = Crc32Lookup[index];
    return ~crc;
}

CRC VS. Fast CRC 過程以 POLY = 1101 舉例如下

                  CRC               Fast CRC 
                                  (用不到 POLY)
====================================================
init CRC        00000000            00000000       
data           ^11010011           ^11010011         
                --------            --------
update CRC      11010011            11010011  -> index
shift           01101001                .
POLY           ^    1101                .
                --------                .
update CRC      01100100                .
shift           00110010                .
shift           00011001                .
                                        .
shift           00001100                .
POLY           ^    1101                .
                --------                .
update CRC      00000001                .
shift           00000000                .
POLY           ^    1101                .
                --------                .
update CRC      00001101                .
shift           00000110                . 
POLY           ^    1101                .
                --------                .
update CRC      00001011                .
shift           00000101                .
POLY           ^    1101                .
                --------                .
update CRC      00001000                .
shift           00000100            00000100  -> lut[index]

Half-Byte Algorithm

參考 Half-Byte Algorithm

雖然可以透過查表加速運算過程,但是這 256 個 32 位元的數值需要龐大的記憶體空間,所以這邊又進化到使用 Half-Byt Algorithm 節省使用記憶體空間

uint32_t crc32_halfbyte(const void* data, size_t length, uint32_t previousCrc32 = 0)
{
  uint32_t crc = ~previousCrc32;
  unsigned char* current = (unsigned char*) data;
  static uint32_t int lut[16] =
  {
    0x00000000,0x1DB71064,0x3B6E20C8,0x26D930AC,0x76DC4190,0x6B6B51F4,0x4DB26158,0x5005713C,
    0xEDB88320,0xF00F9344,0xD6D6A3E8,0xCB61B38C,0x9B64C2B0,0x86D3D2D4,0xA00AE278,0xBDBDF21C
  };
  while (length--)
  {
    crc = lut[(crc ^  *current      ) & 0x0F] ^ (crc >> 4);
    crc = lut[(crc ^ (*current >> 4)) & 0x0F] ^ (crc >> 4);
    current++;
  }
  return ~crc;
}

一次要看 8 個位元太耗記憶體,一次只看 4 個位元的查表就只有 16 個值,但是速度太慢,所以 Half-Byte 演算法採用了每個 byte 拆成兩個 4 個位元處理的設計,每次對一個 byte 的低 4 個位元和高 4 個位元各查表一次,這樣就和一次看 8 位元是一樣的了,對照以下示意圖可以看到 Half-byte Algorithm 兩次查表就對應傳統 CRC 在處理到第 4 位元和第 8 個位元時

  • 因為本來 data 就是一直往右移,所以可以做到一半再次查表
                  CRC           Fast CRC              Half-byte Algorithm            i
=======================================================================================
init CRC        00000000        00000000                  00000000
data           ^11010011       ^11010011                 ^    0011  data_low
                --------        --------                  --------
update CRC      11010011        11010011  index           00000011  crc_low           
shift           01101001           .                         .
POLY           ^    1101           .                         .
                --------           .                         .
update CRC      01100100           .                         .                       0
shift           00110010           .                         .                       1
shift           00011001           .                         .                       2
                                   .                         🡓      lut[crc_low]
shift           00001100           .                      00001100                    
POLY           ^    1101           .                     ^    1101  data_high
                --------           .                      --------
update CRC      00000001           .                      00000001  crc_high         3 
shift           00000000           .                         .
POLY           ^    1101           .                         .
                --------           .                         .
update CRC      00001101           .                         .                       4
shift           00000110           .                         .                       
POLY           ^    1101           .                         .
                --------           .                         .
update CRC      00001011           .                         .                       5
shift           00000101           .                         .
POLY           ^    1101           .                         .
                --------           .                         .
update CRC      00001000           🡓     lut[index]          🡓     lut[crc_high]     6 
shift           00000100        00000100                  00000100                   7

因為一次只看 4 個位元,查表總共有

24=16 個值,和原本 256 個差了 16 倍,和原本建立方式一樣,只是每個值是處理 4 個位元的結果

void generate()
{
    for (size_t i = 0; i < 16; i++) {
        uint32_t crc = i;
        for (uint32_t j = 0; j < 4; j++)
            crc = (crc >> 1) ^ (-(int)(crc & 1) & POLY);
        printf("    0x%08X%s", crc, (i % 4 == 3) ? ",\n" : ", ");
    }
}

所以想要用 Half-byte Algorithm 來實作 CRC-32C 的建表方式只要修改多項式 POLY0x82f63b78,並進行兩次查表即可

Linux 核心中 CRC32 的應用

CRC 案例

挑選 Linux 核心使用 CRC 的案例來說明 (執行 grep -r " CRC " *| less 命令),可說明以下 (用原始程式碼搭配)

輸入grep -r " CRC " *| less 命令後可以看到在 Linux 中尤其在 driver 有大量使用 CRC 的地方

網路七層

image
(圖片來源)
CRC 就是落在 Data Link Layer確保實體傳輸的完整性

e1000e

是 Intel 有線網卡的驅動程式(driver) ,就是在處理透過 PCIe 接上來的 Intel 網卡。當中的netdev.c程式碼就是在處理:

  • 封包的接收與送出(Rx / Tx)
  • 檢查有沒有錯誤(CRC的其他種)
  • 配置網卡參數

e1000e driver

Ethernet 協定有規定 CRC 的除數就是用以下多項式:

P(x)=x³²+x²+x²³+x²²+x¹+x¹²+x¹¹+x¹+x+x+x+x+x²+x+1

因此傳遞的過程是

  • 傳送端會產生 CRC 附在封包的最後 4 個位元
  • 接收端的網卡硬體會驗證 CRC 是否正確
    • 正確:送給 driver 去做其他驗證
    • 不正確會直接把封包丟掉,但會紀錄在CRCERRS後續由 driver 來讀取來統計錯誤
  • driver 會定期讀取硬體錯誤統計資料,了解連線品質,做診斷或告警用

首先e1000_configure_rx 會根據封包大小來決定後續如何處理,有三種分別有三種

if (adapter->rx_ps_pages) {
    /* this is a 32 byte descriptor */
    rdlen = rx_ring->count *
        sizeof(union e1000_rx_desc_packet_split);
    adapter->clean_rx = e1000_clean_rx_irq_ps;
    adapter->alloc_rx_buf = e1000_alloc_rx_buffers_ps;
} else if (adapter->netdev->mtu > ETH_FRAME_LEN + ETH_FCS_LEN) {
    rdlen = rx_ring->count * sizeof(union e1000_rx_desc_extended);
    adapter->clean_rx = e1000_clean_jumbo_rx_irq;
    adapter->alloc_rx_buf = e1000_alloc_jumbo_rx_buffers;
} else {
    rdlen = rx_ring->count * sizeof(union e1000_rx_desc_extended);
    adapter->clean_rx = e1000_clean_rx_irq;
    adapter->alloc_rx_buf = e1000_alloc_rx_buffers;
}
skb 示意圖

image

最常見的路徑會走e1000_clean_rx_irq() ,會先檢查 CRC 最後四個檢查碼是否有被刪除

/* adjust length to remove Ethernet CRC */
if (!(adapter->flags2 & FLAG2_CRC_STRIPPING)) {
    /* If configured to store CRC, don't subtract FCS,
     * but keep the FCS bytes out of the total_rx_bytes
     * counter
     */
    if (netdev->features & NETIF_F_RXFCS)
        total_rx_bytes -= 4;
    else
        length -= 4;
}

CRC 中的平方根運算

開平方根 (square root, sqrt) 的整數運算、原理、加速機制,以及挑選 Linux 核心使用 CRC 的案例來說明,例如:

透過 Linux 核心實例理解 CRC-CCITT(以 rt2800 驅動為例)
rt2800lib.c 是無線網卡的驅動程式,支援 Ralink RT2800 系列晶片

當中的 rt2800_check_firmware_crc(const u8 *data, const size_t len) 函式有用到crc_ccitt 驗證資料

P(x)=x16+x12+x5+1

首先在下列程式碼可以發現,最後兩個 bytes 為 CRC 驗證碼,所以先取出放入fw_crc,因為該資料傳遞方式會 Big-Endian
參考 如何判斷Big-Endian or Little-Endian

static bool rt2800_check_firmware_crc(const u8 *data, const size_t len)
{
	u16 fw_crc;
	u16 crc;
    
	fw_crc = (data[len - 2] << 8 | data[len - 1]);

	crc = crc_ccitt(~0, data, len - 2);

	crc = swab16(crc);

	return fw_crc == crc;
}

然後參考 linux/lib/crc-ccitt.c 後,可以發現實際執行 crc 的是 crc_ccitt_byte
crc_ccitt(u16 crc, u8 const *buffer, size_t len)

u16 crc_ccitt(u16 crc, u8 const *buffer, size_t len)
{
	while (len--)
		crc = crc_ccitt_byte(crc, *buffer++);
	return crc;
}
EXPORT_SYMBOL(crc_ccitt);

再接著參考 linux-xlnx/include/linux/crc-ccitt.h中的crc_ccitt_byte(u16 crc, const u8 c)

static inline u16 crc_ccitt_byte(u16 crc, const u8 c)
{
	return (crc >> 8) ^ crc_ccitt_table[(crc ^ c) & 0xff];
}

可以發現在crc_ccitt_table[(crc ^ c) & 0xff]是利用查表的方式計算 CRC,在 CRC-CCITT 裡,使用的

POLY0x1021,為 16-bit,但一次處理是 8 個位元,所以查表的索引會是 8 個位元,然後回傳的 CRC 會是 16 個位元,如下放上一小段crc_ccitt_table 的內容

u16 const crc_ccitt_table[256] = {
	0x0000, 0x1189, 0x2312, 0x329b, 0x4624, 0x57ad, 0x6536, 0x74bf,
	0x8c48, 0x9dc1, 0xaf5a, 0xbed3, 0xca6c, 0xdbe5, 0xe97e, 0xf8f7, ...}

所以最後再把crc = swab16(crc)進行反轉再和原本的fw_crc做比對
最後使用 swab16(crc) 把 16-bit 的 CRC 驗證碼進行 Byte 交換也就是 Endian 轉換,因為 Firmware 裡儲存的 CRC 與計算出來的 CRC Endian 順序不同

高位和低位交換,例如:
0x12 34 → 0x34 12

另外在以下的核心程式碼中有找到使用CRC的部份

safe_serial.c

drivers/usb/serial/safe_serial.c
在每個 USB 傳輸封包的末尾加上一段包含資料長度與 CRC10 驗證碼
一開始 程式碼 先定義 CRC

#define CRC10_INITFCS  0x000
#define CRC10_GOODFCS  0x000
#define CRC10_FCS(fcs, c) ((((fcs) << 8) & 0x3ff) ^ crc10_table[((fcs) >> 2) & 0xff] ^ (c))

static inline __u16 fcs_compute10(unsigned char *sp, int len, __u16 fcs)
{
	for (; len-- > 0; fcs = CRC10_FCS(fcs, *sp++));
	return fcs;
}
  • fcs 是目前的 CRC 狀態。
  • c:新的 byte
  • (fcs << 8) & 0x3ff:左移 8 bit 並保留 10 bits
  • crc10_table[((fcs >> 2) & 0xff)]:查表,index 是 fcs 的中間 8 bits
  • ^ (c):XOR 上新的輸入 byte
static const __u16 crc10_table[256] = {
	0x000, 0x233, 0x255, 0x066, 0x299, 0x0aa, 0x0cc, 0x2ff,
	0x301, 0x132, 0x154, 0x367, 0x198, 0x3ab, 0x3cd, 0x1fe, ...}

static int safe_prepare_write_buffer(struct usb_serial_port *port, void *dest, size_t size) 中,是 CRC 附加的過程

/* set count */
buf[pkt_len - 2] = count << 2;
buf[pkt_len - 1] = 0;

/* compute fcs and insert into trailer */
fcs = fcs_compute10(buf, pkt_len, CRC10_INITFCS);
buf[pkt_len - 2] |= fcs >> 8;
buf[pkt_len - 1] |= fcs & 0xff;

要放 10-bit CRC,所以要把它:

  • 放在最後 2 個 byte 裡(16 bits),拆成:
    高 2 bits → 放到倒數第 2 byte 的低位
    低 8 bits → 放到倒數第 1 byte
步驟 說明
count << 2 預留最後兩個 bytes,先放 count(協議需求)
fcs_compute10() 用 CRC10 算整包的檢查碼
fcs >> 8 高 2 bit,放進倒數第 2 byte 的低位
fcs & 0xff 低 8 bit,放進倒數第 1 byte

Linux 核心如何使用硬體加速的 CRC32

CRC 原理、如何加速 CRC 運算

TODO: 第 6 週測驗題

測驗題 1, 2 和延伸問題

log
計算

將程式碼大致分為三個區塊

  • acc = 0 : 所需遞迴深度已經完成,直接返回 0
  • num >= base : 代表 num 還可以繼續除以 base,可以繼續遞迴
  • num < base : 計算num 的冪(指數為 base),並繼續進行遞迴

直到遞迴到達指定深度回傳對數,此程式碼對於輸入的限制是 basenum 都是大於 1 的任一浮點數

double logarithm(double base, double num, int acc)
{
    // 1.
    if (acc <= 0)
        return 0.0;
    // 2.
    if (num >= base)
        return 1.0 + logarithm(base, num/base, acc-1);
    // 3.
    double tmp = 1.0;
    for (int i = 0; i < (int) base; i++)
        tmp *= num;    
    return (1.0 / base) * logarithm(base, tmp, acc-1);
}

改寫為更有效的方式

首先可以將計算

numbase 的部分使用快速冪計算,以下用
313
舉例,可以將指數部分
13
轉為二進位,再用指數的二進位展開來計算冪

13=(1101)2=23+22+20

313=383431=1594323

實作如下,我們將指數 exp 當作二進位來看,跌代每個位元,如果該位元是 1,就可以將 num 對應的冪乘到 result 中,直到 exp 為 0

double fast_pow(double num, int exp) {
    double result = 1.0;
    while (exp > 0) {
        if (exp & 1)
            result *= num;
        num *= num;
        exp >>= 1;
    }
    return result;
}

如此可以將原本計算

numbase 時間複雜度
O(base)
減少到
O(log base)
,整體計算對數的時間則是從
O(accbase)
下降至
O(acclog base)

接下來為了讓任意正浮點數皆可以計算,這時候我們又想到了泰勒展開式,並且參考 fdlibm/e_log.c

關於

ln(x) 計算,將
x
拆解成以下

x=2k(1+f)

TODO:使用定點數改寫

PELT 中的指數函數

指數函數的計算 (特別是定點數)、原理,挑選 Linux 核心使用指數函數的案例,如 PELT 用到的 EWMA

PELT(Per-Entity Load Tracking)

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.

由於是每個週期是

1024 毫秒,並且需要一個衰減因子
y
經過
32
次週期後剩一半

y32=0.5

這個部分就是使用指數衰減的展現,以下說明指數加權移動平均 (EWMA)

Exponentially Weighted Moving Average (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

SMAk=1ki=nk+1npi

EWMA 會讓資料隨著時間權重降低,其中 smoothing factor 平滑因子

α 決定這個資料每次要衰減多少,且
0α1

exponential functions are used to assign exponentially decreasing weights over time

s0=x0st=αxt+(1α)st1,t>0

EWMA 會讓最新的資料

xt 有較高的權重
α
,其餘過去資料會因為不斷乘上
(1α)
而變得越來越小,如下
s0=x0s1=αx1+(1α)s0=αx1+(1α)x0s2=αx2+(1α)s1=αx2+(1α)αx1+(1α)2x0s3=αx3+(1α)s2=αx3+(1α)2αx1+(1α)3x0st=αxt+(1α)st1=αxt+(1α)αxt1+(1α)2αxt2++(1α)t1αx1+(1α)tx0=α[xt+(1α)xt1+(1α)2xt2+(1α)3xt3++(1α)t1x1]+(1α)tx0

隨著時間

t 可以忽略
x0
這一項

st=αk=0(1α)kxtk

這邊的

(1α) 即可視為 PELT 所用到的衰減因子
y
x
的部分用
u
代替方便對應程式碼說明

st=αk=0ykutk

Linux Kernel Document: Schedutil 文檔所提及的 EWMA 的計算將與幾何級數成正比的部分獨立出來成為一個函數

ewma_sum(u)

ewma_sum(u)=k=0ykutk

完整的 EWMA 就可以是

ewma_sum(u) 和 幾何級數
ewma_sum(1)
的比值,這樣也可以減少原本因為無窮級數的計算丟失精度

ewma(u)=ewma_sum(u)ewma_sum(1)

將分母與分子展開,可以得到我們一開始推導 EWMA 的公式,因為

|y|<1 ,所以分母
ewma_sum(1)
部分即為幾何級數,最後一定會收斂所以可以寫成
11y

ewma(u)=k=0ykutkk=0yk=k=0ykutk11y=k=0ykutk1α=αk=0ykutk

程式碼與其對應數學意義

PELT 原始碼參考: linux/kernel/sched/pelt.cpelt.hsched-pelt.hsched.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.hsched_avg 用於定義一個計算負載平均的物件,其中跟 running 有關的數值就是 util_sum (運用CPU的時間),與 runnable 有關的數值是 runnable_load_sum (等待的時間),而負載是當前時間任務對系統產生的壓力,也就是 load_sum 應是 running + runnable 的時間
TODO: 這裡要計算的各種指標代表意義

struct sched_avg {
	u64				last_update_time;
	u64				load_sum;
	u64				runnable_sum;
	u32				util_sum;
	u32				period_contrib;
	unsigned long			load_avg;
	unsigned long			runnable_avg;
	unsigned long			util_avg;
	unsigned int			util_est;
} ____cacheline_aligned;

decay_load

目的為計算

val×yn

val 單位是?

也就是每個值每次衰減的計算,

n 代表第幾次衰減,我們期望最終在 負載平均週期數 LOAD_AVG_PERIOD
32
次衰減後剩原本的
12

y32=12

因為 kernel 中沒有 FPU , 直接算

yn 較複雜,所以將
yn
整合成以
12
為底數的冪,就可以用右移的位元操作完成這部分的指數運算

yn=(12)n32

為了避免

n32 因無法整除而丟失精度,首先觀察
n32
,商
k
為整數

n=k32+n%32

利用指數性質

yn=y(a+b)=ya×yb

yn=(12)k3232×(12)n%32=(12)k×(12)n%32

這樣就可以將

val 右移
k
,並用查表方式來取得
(12)n%32
,實作可以看到 pelt.c#L49

if (unlikely(local_n >= LOAD_AVG_PERIOD)) {
    val >>= local_n / LOAD_AVG_PERIOD;
    local_n %= LOAD_AVG_PERIOD;
}
val = mul_u64_u32_shr(val, runnable_avg_yN_inv[local_n], 32);

其中 runnable_avg_yN_inv 即為查詢所有可能

(12)n%32 的查表(將浮點數轉成 32 bit 的二進位結果),再用 mul_u64_u32_shr
val×(12)k
(12)n%32
這兩個值相乘取得最後衰減後的數值

include/linux/math64.h 有支持 128 位元整數和不支持的 mul_u64_u32_shr 兩種計算方式,如果沒有支持 128 位元整數的編譯器,可以用定點數的方式計算,實作如下

static __always_inline u64 mul_u64_u32_shr(u64 a, u32 mul, unsigned int shift)
{
    u32 ah = a >> 32, al = a;
    u64 ret;

    ret = mul_u32_u32(al, mul) >> shift;
    if (ah)
        ret += mul_u32_u32(ah, mul) << (32 - shift);
    return ret;
}

用定點數來解釋這段乘法,首先 a 拆解成高 32 位元的 ah 和 低 32 位元的 al

a=ah232+al

amul 的相乘即可分成高位元和低位元各自乘上 mul 後的和

amul=(ah232+al)mul=ahmul232+almul

先計算

almul

ret = mul_u32_u32(al, mul);

再加上

ahmul232

ret += mul_u32_u32(ah, mul) << 32;

shift 在這裡作用是甚麼 ? 如果看 shift 在整個衰減計算的位置,相當於

yn=(12)k×(12)n%32×(12)shift

也就是可以防止衰減太快,有縮放衰減速度的作用

accumulate_sum

將所需計算的資料分為三個階段,d1 是上一個週期中剩下的時間,d1 是有經過完整週期的時間段的數量,d3 是目前週期已經過的時間,參考 Yiwei Lin 的 Linux 核心設計: Scheduler(4): PELT 說明 accumulate_sum

           d1          d2           d3
           ^           ^            ^
           |           |            |
         |<->|<----------------->|<--->|
 ... |---x---|------| ... |------|-----x (now)

所以計算這三個階段總和

ewma_sum(u)
u=uyp+d1yp+1024n=1p1yn+d31

加總所有時間到

delta ,包括過去時間片段 sa->period_contrib

delta += sa->period_contrib;

其中

p 是目前已經過完的週期數量,一個週期是
1024
微秒,如果已經過了
delta
微秒,則
p

p=delta1024

Step 1. 計算過去累積的負載

u 目前的衰減後的數值,即
uyp
,程式碼如下,用 decay_load 更新累積至今的總和 *_sum

sa->load_sum = decay_load(sa->load_sum, periods);
sa->runnable_sum = decay_load(sa->runnable_sum, periods);
sa->util_sum = decay_load((u64)(sa->util_sum), periods);

Step 2. 先更新 delta 成無法組成一個周期的剩餘時間,如果當前沒有負載的話可以不用計算

d1yp+1024n=1p1yn+d31 ,如果有則交由 _accumulate_pelt_segments 計算

/*
 * Step 2
 */
delta %= 1024;
if (load) {
    /*
     * This relies on the:
     *
     * if (!load)
     *	runnable = running = 0;
     *
     * clause from ___update_load_sum(); this results in
     * the below usage of @contrib to disappear entirely,
     * so no point in calculating it.
     */
    contrib = __accumulate_pelt_segments(periods,
            1024 - sa->period_contrib, delta);
}

最後加上新的負載貢獻,根據

ewma_sum(u),最新的貢獻值不用乘上衰減因子,直接加上去就好

if (load)
    sa->load_sum += load * contrib;
if (runnable)
    sa->runnable_sum += runnable * contrib << SCHED_CAPACITY_SHIFT;
if (running)
    sa->util_sum += contrib << SCHED_CAPACITY_SHIFT;

__accumulate_pelt_segments

計算當前負載貢獻值,即 d1d2d3 三個部分 ,也就是

d1ypc1+1024n=1p1ync2+d31c3

d1 部分,即

d1yp ,有多少貢獻值就直接乘上衰減因子

c1 = decay_load((u64)d1, periods);

d2 部分,即

1024n=1p1yn,由於這個計算會需要迴圈,這裡進行拆解以方便計算,

n=1p1yn=n=0yn第一部分n=pyn第二部分y0第三部分

第一部分

n=0yn 就是無窮幾何級數,因為
|y|<1
最後會收斂

n=0yn=11y

1024×n=0yn 會是一個常數,所以定義此數 LOAD_AVG_MAXlinux/kernel/sched/sched-pelt.h,估算
y
如下

y=0.51320.97857

所以

LOAD_AVG_MAX 大約是
47788

LOAD_AVG_MAX=1024×11y=1024×110.9785747788

有一些誤差 LOAD_AVG_MAX 是定義 47742

第二部分

n=pyn 展開如下

n=pyn=yp+yp+1+yp+2+=yp×(1+y+y2+)=yp×n=0yn

乘上

u=1024 的結果,即使用 decay_load 來計算
yp×LOAD_AVG_MAX

1024×yp×n=pyn=yp×(1024n=0yn)=yp×LOAD_AVG_MAX

第三部分即

1024×y0=1024

d2 的計算完整如下

c2=1024n=1p1yn=1024[n=0ynn=pyny0]=1024[n=0yn(ypn=0yn)1]=1024[(11y11yyp)1]=10241y10241yyp1024=LOAD_AVG_MAXLOAD_AVG_MAXyp1024

所以 d2 實作如下

c2 = LOAD_AVG_MAX - decay_load(LOAD_AVG_MAX, periods) - 1024;

d3 部分,最新進的貢獻,不用乘上衰減因子

最終回傳當前負載

return c1 + c2 + d3;

___update_load_sum

___update_load_avg

以上函數都是和

ewma_sum(u) 相關的計算,完整的
ewma(u)
最後還需除以
ewma_sum(1)

ewma(u)=ewma_sum(u)ewma_sum(1)

首先呼叫 get_pelt_divider() , 相當於取得

ewma_sum(1)

___update_load_avg(struct sched_avg *sa, unsigned long load)
{
	u32 divider = get_pelt_divider(sa);
...

最後相除,即完成 EWMA 的計算

sa->load_avg = div_u64(load * sa->load_sum, divider);
sa->runnable_avg = div_u64(sa->runnable_sum, divider);

get_pelt_divider

依照原先推倒的公式來說,應該是直接回傳

ewma_sum(1)=n=0yn

#define PELT_MIN_DIVIDER	(LOAD_AVG_MAX - 1024)

static inline u32 get_pelt_divider(struct sched_avg *avg)
{
	return PELT_MIN_DIVIDER + avg->period_contrib;
}

誤差分析

Linux 核心中的浮點數運算

2 的冪乘法

實作 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 my_ldexpf(float x, int exp);

初步想法

先將 float 轉成 unsigned 然後再去取 signexponentmantissa

float my_ldexpf(float x, int n) // x * 2^{n} ==相當於: x <<= n
{		
    int sign = (unsigned)x >>31;
    int exponent = ((unsigned)x>>23)& 0xff;
    int mantissa = (unsigned)x & 0x7fffff;
}

根據 IEEE-754 Floating Point Converter 表示法去嘗試當 x 是 -1.2345 印出 signexponentmantissa 但除了 sign 其他都不對,如下

     sign  exponent   mantissa
正確: 1      127       1967129
我的: 1      255       2097151   

檢討
關於浮點數和整數的轉換可以看到以下 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 就會只剩整數位的

1,雖然可以取到正確的 sign ,但根本取不到正確的 exponentmantissa

 expression   decimal   binary      
          x   -1.2345   10111111 10011110 00000110 00010000
(unsigned)x   -1        11111111 11111111 11111111 11111111

改進本筆記的書寫,讓排版接近指定教材的風格

分析

了解浮點數和整數轉換的問題後,我們知道不行直接將浮點數用 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] 有解釋到:

問題

  • 為什麼這跟 type punning 有關?
  • 「由沒有字元類型的左值表達式讀取」甚麼意思?
  • 「如果這種表示是由副作用產生的」是指甚麼副作用?
    side effect: 有 an unexpected result of a situation 未曾預料的結果 的意思

出自 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 浮點數後輸出的

11111111 11111111 11111111 11111111 是代表
1
還是未定義 NaN

TODO: 所以是怎麼取值的

所以將浮點數轉為二進位數的方式如下

union{
    float f;
    uint32_t u32;
}val;

val.f = x;

取出 sign, exp, frac;

uint32_t  sign = val.u32 >> 31;
uint32_t  exp = (val.u32 >> 23) & 0xff;
uint32_t  frac = val.u32 & 0x7fffff;

轉換為二進位置後,就可以來看要怎麼進行相乘運算

浮點數的二進位制運算

更多浮點數的運算可以看到教材 CS:APP 第 2 章重點提示和練習,所以取得 signexpfrac 就可以進行單精度浮點數和

2n 的乘法運算了,以下使用
x=1.2345
為例

首先浮點數編碼成二進位的方式為

value=(1)s×M×2E

也就是我們只要找到

s (sign) 、
M
E
就可以運算浮點數和
2n
的乘法,
s
已經知道是
1
,將
1.2345
轉換成二進位制會是無限循環,由於是單精度浮點數,所以就取到
23
位即可

1.2345=(1)1×1.00111100000010000011001...2×20

  • Significand:

    M=1.001111000000100000110012frac=00111100000010000011001

  • Exponent

    E=0Bias=127exp=127=01111111

  • Result

    
    
    
    
    
    
    binaryGraph
    
    
    
    sign
    
    sign
    
    1
    
    
    
    exp
    
    exp
    
    01111111
    
    
    
    
    frac
    
    frac
    
    00111100000010000011001
    
    
    
    
    

現在回到題目,我們這個時候會取得需要相乘的 2 的

n 次方,而 2 的次方 2 的冪 之間的相乘就是將次方數直接相加,如下
2E×2n=2E+n

而此時的

exp 會變成
exp=E+Bias

exp=(E+n)+Bias

那其實剛剛已經取出 signexpfrac,所以可以直接將 n 加到 exp

int new_exp = exp + n;

那最後的答案就會是這三個值的 logical or

val.u = sign << 31 | new_exp << 23 | frac

完全不會動到 fracsign,只要處理 exp 的部份,其實只要理解浮點數在電腦中的編碼模式就能解這題了

最後回傳 val.f 就是答案

測試

CS:APP 第二章 有提到浮點數編碼後會有三種狀態

  • denormalized :
    exp=00...00
  • normalized :
    exp00...00
    and
    exp11...11
  • special :
    exp=11...11

在運算 new_exp 之前可以先檢查 x 是否為有效值

  1. xNaN,其實就是在檢查 exp 的部份

    ​​​​ if (val.u32 > 0x7f800000)...
    
  2. x 是無限大或無限小

    ​​​​if(val.u32 == 0x7f800000)...
    
  3. x 是 0

    ​​​​if(val.u32 == 0)...
    

再來需要檢查 new_exp 是否為有效值

  1. new_exp overflow
    因為
    exp
    只有 8 個位元,所以需要檢查 new_exp 是否超過
    128
    也就是二進位制的
    11111111
    ​​​​if (new_exp >= 0xff)...
    
  2. new_exp underflow
    ​​​​if (new_exp <= 0)...
    
  3. x 是 normal 但輸出為 denormal
    ​​​​if (new_exp < 1)...
    

完整程式碼

使用 lab0-c 規範的程式碼風格 (類似 Linux 核心原始程式碼)。

float my_ldexpf(float x, int n)
{
    union{
        float f;
        uint32_t u32;
    }val;
    
    val.f = x;
    uint32_t sign = val.u32 & 0x80000000;
    uint32_t exp = (val.u32 >> 23) & 0xff;
    uint32_t frac = val.u32 & 0x7fffff;
    
    // x = +-inf or 0 or n=0
    if (exp >= 0xff || x == 0 || n ==0)
        return x;
        
    uint32_t new_exp = exp + n;
    
    if(new_exp >= 0xff){ // overflow or nan
        val.u32 = sign | 0x7f800000;
        return val.f;
    }
    
    if(new_exp <=0){ // underflow
        val.u32 = sign;
        return val.f;
    }
    
    val.u32 = sign | new_exp << 23 | frac;
    
    return val.f;
}

參考: ReactOS 中的 ldexpf 實作

當你列出超連結時,應當給予對應的標題

驗證不夠,且需要更精簡

image
只需檢查 exp 是否發生 overflow
或是 underflow
,當檢查到一位時就將
frac
設置為 0,因為 任何數 與 正負無窮大
做運算皆會得到 NaN (也就是未定義的數學結果)

#include <stdint.h>
float my_ldexpf(float x, int n)
{
    union{
        float f;
        uint32_t u32;
    }val;
    
    val.f = x;
    uint32_t sign = val.u32 & 0x80000000;
    uint32_t exp = (val.u32 >> 23) & 0xff;
    uint32_t frac = val.u32 & 0x7fffff;
            
    exp += n;
    if(exp < 0 || exp >= 0xff){        
        frac = 0; 
        exp = 0xff;
    }           
    
    val.u32 = sign | exp << 23 | frac;
    
    return val.f;
}

幾何平均數

要求:實作 geomean,計算幾何平均數,不能用 FPU

int geomean(float *in, int n) { } // n 表示 in 的總數

GM=val1val2valnn

補足背景知識

之前已探討單精度浮點數的格式,此處可略過。

首先根據 IEEE 754 單精度浮點數佔用 32 位元的空間,其編碼方式為:

value=(1)s×M×2E
image

  • fraction(M)部分存放小數點後資料
  • 舉例: 8.5 = 2³ + 1/2¹ 轉成二進制為 1000.1 調整成 1.0001 x 2³
    • E的部分就是3,但exp編碼以01111111 (127)表示指數為0,所以指數為3時應先+3,表示成10000010 (130)
    • frac(M) 的部分就是0001因此放到 mantissa 裡就是00010000000000000000000 (共23bits)
  • uint32_t a的三個部分取出方法
    ​​​​uint32_t sign_a = a >> 31;
    ​​​​uint32_t exp_a  = (a >> 23) & 0xFF; //0xFF = 1111 1111₂ 把以外的位元都清0
    ​​​​uint32_t frac_a = a & 0x7FFFFF;
    

將浮點數轉成定點數時,涉及數值誤差,應予以討論。

TOOD: 數值誤差

定點數

維基百科 介紹中描述了定點數表示法中,整數和小數的部分會用相同底數

b 來表示(就是會是一樣是十進位或是二進位),如果此有禮數 有理數的小數的部分到第
n
位,此有理數會是
bn
的整數倍

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

bn.

乍看之下,其實跟浮點數很相似,都是表示有理數的一種方式,但是可以看到,定點數會有 scaling factors 可讓數值表示成整數在運算,如下:

假設我們現在要將

0.123
2.5
相乘,我們可以想成這兩個定點數相乘
(123×103)×(25×101)

縮放因子在此次運算中即為

103×101=104

此時就可以不用用到 FPU 的將兩個整數相乘

123×25=3075

最後再去乘上縮放因子即可得到結果

3075×104=0.3075

回到 Linux 核心中的浮點數相乘,我們如果可以用定點數的概念,那麼就可以用以下步驟得到結果

  1. M
    變成定點數
  2. M
    整數部分相乘
  3. E
    和縮放因子的指數做簡單的加法運算

浮點數相乘

接著兩個單精度浮點數相乘

a×b
value=(1)sa+sb×(Ma×Mb)×2Ea+Eb

  • sign的部分其實就是把
    sa
    sb
    做XOR
  • 接下來
    M
    的部分因為一開始在存放時,已省略 1.0001 最高位元的 1 ,所以在做乘法時應當要補回 mantissa
    M
    最前面 最高位元的 1: M_a = (1 << 23) | frac_a; 再做相乘
  • exp 的部分就是直接相加

注意用詞的精準,避免說「最前面」,而是「最高位元」

除了以上大概念外,還會遇到像是 mantissa 相乘時會有 overflow 和正規化的問題的問題

  • IEEE‑754 要求 mantissa 必須是 1.xxxx 的形式,也就是在 1.0~2.0 之間。如果乘出來 ≥ 2.0,就要把尾數右移一位,並把 exponent 加 1。

  • 參考你所不知道的 C 語言: 浮點數運算中有提到浮點數的加法,若超過1也是把其中一個 mantissa 往右移,然後 exp 加一。

    image

  • 所以乘法應當也會有兩種情況

    乘積大小 Bit47 代表值範圍 處理方式
    ≥ 2.0 1 [2.0, 4.0) P >>= 24; exp++
    [1.0,2.0) 0 [1.0, 2.0) P >>= 23;
    ​​​​uint64_t P = M_a * M_b;
    
    ​​​​// 正規化:若 ≥2.0,右移 24, exp+1;否則右移 23
    ​​​​if (P & (1ULL<<47)) { //假設Bit47是1則代表乘積超過2
    ​​​​    P >>= 24;
    ​​​​    exp_r += 1;
    ​​​​} else {
    ​​​​    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 的情況

floating point precise

參考: 留意浮點數運算的陷阱to nearest, ties to even作法

第一個方法:
M
先相乘再轉換成 2 的冪

M 全部相乘,並使用
log
取指數再進行開根號運算,如下:

初始化一個

Mtotal ,假設
Mtotal
有 64 個位元
Mtotal=1

對傳入的陣列

in 進行迭代取出其數值的
fraci
部分,則其
Mi
的定點數為
Mi=(223+fraci)×223

Mtotal 會相乘每個
Mi
轉成定點數的數值部分
Mtotal=MtotalMival

最後開根號

GM=[(Mtotal223n)×2Esum]1n=[2log2Mtotal223n2Esum]1n=21n(log2Mtotal23n+Esum)

假設陣列

in={1.2,1.3,9.4}

i=0
M0=1001100110011001100110102223

Mtotal=1M0val=1001100110011001100110102

此時的 scaling factors

sf=223

i=1
M1=1010011001100110011001102223

Mtotal=10011001100110011001101021010011001100110011001102=110001101011001101001110100011001011001001001002

為了避免下一次和另一個 24 位元的數值相乘完溢位,需要將 47 位元進行捨入到 40 位元,這裡的捨入採用

Mtotal=11000110101100110100111010001100101100102×27

此時的 scaling factors

sf=22322327=239

i=2
M2=1001011001100110011001102220

Mtotal=110001101011001101001110100011001011001021001011001100110011001102=10010100111100110100111000101100101100100000000011001100110011002260

此時的 scaling factors

sf=239220=259

最後開根號,因為浮點數的整數部分都是 1 ,所以

Esum 會是 0
GM=259+60212

似乎結果和第二個方法一樣,在沒有 FPU 的條件下可以只看

E 的部分就好

int round(M)
{    
    int mult_2 = M & (1 << 49);
    
}

int geomean(float *in, int n) 
{
    union{
        float f;
        uint32_t u32;
    }converter;
    
    int M_total = 1, error = 0;
    uint32_t e_sum = 0, sf = 0;
    
    for(auto num: in){
        converter.f = num;
        uint32_t exp = (converter.u32 >> 23) & 0xff;
        uint32_t frac = converter.u32 & 0x7fffff;
        
        uint32_t M_val = (1 << 23) | frac;
        
        M_total *= M_val;
        
        error = round(M_total);
        
        e_sum += exp;        
    }
    
    e_sum +=  - (n * 23)    
    
    
}

第二個方法: 將
M
之間的相乘轉成 2 的冪之間的相乘

原先上面認為

M 要先全部相乘完再去做
log
取指數 ,但下面發現也能先對各別的
M
log
完再去相乘,因此我們想比較哪一個作法會保留更多的精準度。

如此,幾何平均數的運算可以變成 2 的冪的運算,所以在浮點數編碼成二進位後,就只要做指數部分的運算了

GM=[(M0M1Mn)×(2α02α12αn)]1n=[(2β02β12βn)×(2α02α12αn)]1n=21ni=0n(αi+βi)

浮點數的

M 要轉換成底數為 2 的冪需要算出指數部分
β

M=2ββ=log2M

Exponent 的總結果

Etotal 會是
Etotal=1ni=0nαi+βi=1ni=0n(αi+log2Mi)=1ni=0n[αi+log2(frac223+1)]1ni=0nαi

接著會發現

log2(frac223+1) 的部分計算結果恆為零,因為
frac
本來就只有 23 個位元,乘上
223
就相當於右移 23 個位元直接變成 0,所以只剩
1ni=0nαi
需要計算。

運算 $log_2M$ 的部分則可以透過以下實作完成,但是就沒辦法很精準,因為只能算出整數
int log = -1; 
while(M) { 
    log++; 
    M >>= 1; 
} 

參考: What's the best way to calculate log(base2) of any number?

實作

根據剛才的分析, log2 的部分就不用算了

int log2(uint32_t frac)
{
    int log = -1;    
    while(frac) { 
        log++; 
        frac >>= 1; 
    }
    return log < 0 ? 0 : log;
}

int geomean(float *in, int n) 
{ 
    union{
        float f;
        uint32_t u32;
    }val;
    
    uint32_t e_sum = 0;
    
    for(int i=0; i<n; i++){
        val.f = in[i];
        uint32_t exp = (val.u32 >> 23) & 0xff;
        uint32_t frac = val.u32 & 0x7fffff;

        e_sum += exp; // + log2(frac / (1<<23) + 1);
    }
    e_sum /= n;
    
    if(exp_sum < 0 || exp_sum >= 0xff){        
        frac = 0; 
        exp_sum = 0xff;
    } 

    val.u32 = e_sum << 23;
    return (int)val.f;
} 

第一個和第二個方法都沒有辦法到很精確的數值,因為

frac 的部份最後都會被捨棄,即使取
log
也沒有辦法到很精確,所以可以用比較直觀的第二個方法實做

誤差分析

TODO: 如何比較兩種浮點數操作方法的精準度損失?
為了比較兩種實作的精準度,我們使用 double 作為 ground truth,接著比較兩種作法它們和 ground truth 的誤差,所以使用相對誤差作為指標:error = |approx - truth| / |truth|

expm1(x) 討論

延伸 BigMickey expm1 相關探討

此處討論 linux kernel 中浮點數運算對 expm1() 的影響

double expm1(x){} // 回傳 exp(x) - 1

exp(x)

要使用單精度浮點數還是雙精度浮點數進行分析?

在 IEEE 754 中對浮點數進行二進位編碼,所以可以將

ex 轉為與 2 的冪相關的數就能更容易計算,參考 ulysses-libc/src/math/exp.c 如下

我們期望

x
k
倍的
ln2
並且加上誤差
r
,所以將 x 分解成

x=kln2+r

所以

k 應近似於

k=round(x1ln2)

ex 則可以分解成

ex=ekln2+r=ekln2×er=2k×er

r
x
kln2
的誤差,也就是一個正確值
hi
與近似值
lo
的差

r=xkln2=hilo

|r|<=0.5×ln20.34657359028

如此可以專注於需要高精度計算的

er 部分

er 和解
ex
的數值範圍並不一樣,這種拆解方式可以防止直接算
ex
時益位

以往計算

ex 都是直接使用泰勒展開式

er=1+r+r22!+r33!+r44!=1+r+r22+r36+r424+

但是電腦沒有無窮大,所以會丟失精度,在參考的程式碼中定義了一個 special rational 函數

R(r2)

R(r2)=r×(er+1)(er1)=r×2+r+r22+r36+r424+r+r22+r36+r424+=r×2r+1+r2+r26+r324+1+r2+r26+r324+

u=r2+r26+r324+ 代入

R(r2)=r×2r+1+u1+u

由幾何級數的概念可推導出:

  • 第一項
    2r/(1+u)

2r1+u=2r×(1u+u2u3+u4+)

  • 第二項
    1/(1+u)

11+u=1u+u2u3+u4+

  • 第三項
    u/(1+u)

u1+u=u×(1u+u2u3+u4+)

所以

R(r2) 前 5 項是

R(r2)=r×[(2r+u+1)(1u+u2u3+u4+)=2+r26r4360+r615120r8604800

此時發現使用雷米茲演算法的五階多項式就能將誤差控制在

259

R(z)2+P1z+P2z2+P3z3+P4z4+P5z5

,

z=r2,而其中的 5 個係數就是

P1=16=1.66666666666666019037e01

P2=1360=6.61375632143793436117e05

P3=115120=6.61375632143793436117e05

P4=1604800=1.65339022054652515390e06

P5=4.13813679705723846039e08

P5 好像不是 1/30,240,000

R(r2) 得到
er

R(r2)=r×(er+1)(er1)er=r+R(2)R(r2)rer=1+2rR(r2)r

R(z) 代入計算

er=1+2rR(z)r

此時討論如果

r 接近 0:

R(r2)2+

R(r2)r2r+

計算

2r2r 時,分母中較大的
2
和較小的
r
相減可能會失去精度,如:

假設支持 15 位的精準度,

r=0.0001:

20.0001+0.000126=1.999900016

但是因為只能支持 15 位,後面的位數直接被捨棄剩

1.999900016666667

所以整理成

er=1+r+rc(r)2c(r)

c(r)=r(P1r2+P2r4+P3r6+P4r8+P5r10)

因為

c(r) 算起來會很小,所以
2c(r)
相對
2r
會更接近 2,減少精度損失,而且可以讓
1+r
單獨計算

算出

er 後,就可以用 my_ldexpf 算出
ex=2k×er

expm1