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

執行人: kuanyu0712, Hlunlun

TODO: 第 2 週測驗題

測驗題 1, 2 和延伸問題

測驗一整理

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

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

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

• mpi_t structure 示意圖:

先整理各個函式的功用

• mpi_init: 初始化 mpi_t 記憶體空間
• mpi_clear: 在整個運算後釋放 mpi_t 裡 data 所指向的所有記憶體
• 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 bits 的 op 先取出最低位的 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 (該例子不會觸發ｃ!=0的情況)
  • 當第一次時 tmp->data[k] += (r & INTMAX) + c 第 k 個位置就會是相乘完後結果 r 的最低 31 個位元
  • 假設前面還有 5 個位元，則ｒ>>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)
   ┌────
11 │1101  ← 被除數 (13)
   │11↓↓  ← 11 ≥ 11，可以除，商寫1，餘數=0
   │ 00↓  ← 帶下來0，得00 < 11，商寫0
   │  01  ← 帶下來1，得01 < 11，商寫0
   │   1  ← 餘數 (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-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=\mathrm{商}+CRC$$

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

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

data: 11010011
POLY: 1011

可以參考 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
              ----
                11             

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

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

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

反轉位元組順序 Swap Endianness

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

反轉 POLY 成 1101，如果最低位是 1 就進行 XOR，如果是 0 舊只右移，重複 8 次，最後得到餘數 100
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 值共有

$2^8=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) ^ POLY) : (crc >> 1);

等同於

crc = (crc & 1) ? ((crc >> 1) ^ POLY) : (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 = 00001101 舉例如下

                  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 個位元時

                  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 個位元，查表總共有

$2^4=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 的建表方式只要修改多項式 POLY 為 0x82f63b78，並進行兩次查表即可

Linux 核心中 CRC32 的應用

• 比較檔案是否相同

• 模組載入與驗證

• 封包錯誤偵測

設計實驗來分析 Fast CRC32

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

實作小型的模組化合成器

Q notation

TCP 連線

TODO: 第 6 週測驗題

測驗題 1, 2 和延伸問題

$log$ 計算

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

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

double logarithm(double base, double num, int acc)
{
    if (acc <= 0)
        return 0.0;
    
    if (num >= base)
        return 1.0 + logarithm(base, num/base, acc-1);

    temp = num^base    
    return (1.0 / base) * logarithm(base, temp, acc-1);
}

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 然後再去取 sign、exponent、mantissa

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 印出 sign、exponent、mantissa 但除了 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 ，但根本取不到正確的 exponent 和 mantissa

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

$11111111111111111111111111111111$ 是代表

$-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 章重點提示和練習，所以取得 sign、exp、frac 就可以進行單精度浮點數和

$2^n$ 的乘法運算了，以下使用

$x=-1.2345$ 為例

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

$$value=(-1{)}^s\times M\times 2^E$$

也就是我們只要找到

$s$ (sign) 、

$M$、

$E$ 就可以運算浮點數和

$2^n$ 的乘法，

$s$ 已經知道是

$1$，將

$1.2345$ 轉換成二進位制會是無限循環，由於是單精度浮點數，所以就取到

$23$ 位即可

$$-1.2345=(-1{)}^1\times {1.00111100000010000011001...}_2\times 2^0$$

• Significand:
  

  $$\begin{array}{rl}M& =1.{\underset{―}{00111100000010000011001}}_2\\ frac& =00111100000010000011001\end{array}$$

• Exponent
  

  $$\begin{array}{rl}E& =0\\ Bias& =127\\ exp& =127=01111111\end{array}$$

• Result

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

$n$ 次方，而 2 的次方 2 的冪 之間的相乘就是將次方數直接相加，如下

$$2^E\times 2^n=2^{E+n}$$

而此時的

$exp$ 會變成

$$exp=E+Bias$$

$$ex{p}^{\prime }=(E+n)+Bias$$

那其實剛剛已經取出 sign、exp、frac，所以可以直接將 n 加到 exp

int new_exp = exp + n;

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

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

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

最後回傳 val.f 就是答案

測試

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

• denormalized :
  $exp=00...00$
• normalized :
  $exp\ne 00...00$ and
  $exp\ne 11...11$
• special :
  $exp=11...11$

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

1. x 是 NaN，其實就是在檢查 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 實作

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

驗證不夠，且需要更精簡

只需檢查 exp 是否發生 overflow

$\mathrm{\infty }$ 或是 underflow

$-\mathrm{\infty }$，當檢查到一位時就將

$frac$ 設置為 0，因為 任何數 與 正負無窮大

$\mathrm{\infty }$ 做運算皆會得到 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 的總數

補足背景知識

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

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

$$value=(-1{)}^s\times M\times 2^E$$

• 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$ 位，此有理數會是

$b^{-n}$ 的整數倍

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

$b^{-n}$.

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

假設我們現在要將

$0.123$ 和

$2.5$ 相乘，我們可以想成這兩個定點數相乘

$$(1230\times {10}^{-3})\times (25\times {10}^{-1})$$

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

$${10}^{-3}\times {10}^{-1}={10}^{-4}$$

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

$$123\times 25=3075$$

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

$$3075\times {10}^{-4}=0.3075$$

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

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

浮點數相乘

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

$a\times b$

$$value=(-1{)}^{s_a+s_b}\times (M_a\times M_b)\times 2^{E_a+E_b}$$

• sign的部分其實就是把
  $s_a$ 和
  $s_b$ 做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 加一。
  

  

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

  乘積大小	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 的情況

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

第一個方法:

$M$ 先相乘再轉換成 2 的冪

將

$M$ 全部相乘，並使用

$log$ 取指數再進行開根號運算，如下:

初始化一個

$M_{total}$ ，假設

$M_{total}$ 有 64 個位元

$$M_{total}=1$$

對傳入的陣列

$in$ 進行迭代取出其數值的

$fra{c}_i$ 部分，則其

$M_i$ 的定點數為

$$M_i=\underset{―}{(2^{23}+fra{c}_i)}\times 2^{-23}$$

用

$M_{total}$ 會相乘每個

$M_i$ 轉成定點數的數值部分

$$M_{total}=M_{total}\cdot M_{i_{val}}$$

最後開根號

$$\begin{array}{rl}GM& =[(M_{total}\cdot 2^{-23n})\times 2^{E_{sum}}‘‘{]}^{\frac{1}{n}}\\ & =[2^{lo{g}_2{M}_{total}}\cdot 2^{-23n}\cdot 2^{E_{sum}}{]}^{\frac{1}{n}}\\ & =2^{\frac{1}{n}(lo{g}_2{M}_{total}-23n+E_{sum})}\end{array}$$

假設陣列

$in=\{1.2,1.3,9.4\}$

當

$i=0$

$$M_0=\underset{―}{100110011001100110011010}2\cdot 2^{-23}$$

$$\begin{array}{rl}{M}_{total}& =1\cdot M_{0_{val}}\\ & ={100110011001100110011010}_2\end{array}$$

此時的 scaling factors

$$sf=2^{-23}$$

當

$i=1$

$$M_1={\underset{―}{101001100110011001100110}}_2\cdot 2^{-23}$$

$$\begin{array}{rl}{M}_{total}& ={100110011001100110011010}_2\cdot {101001100110011001100110}_2\\ & =110001101011001101001110100011001011001001001002\end{array}$$

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

$$M_{total}={1100011010110011010011101000110010110010}_2\times 2^7$$

此時的 scaling factors

$$sf=2^{-23}\cdot 2^{-23}\cdot 2^7=2^{-39}$$

當

$i=2$

$$M_2=\underset{―}{100101100110011001100110}2\cdot 2^{-20}$$

$$\begin{array}{rl}{M}_{total}& ={1100011010110011010011101000110010110010}_2\cdot {100101100110011001100110}_2\\ & ={1001010011110011010011100010110010110010000000001100110011001100}_2\\ & \approx 2^{60}\end{array}$$

此時的 scaling factors

$$sf=2^{-39}\cdot 2^{-20}=2^{-59}$$

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

$E_{sum}$ 會是 0

$$\begin{array}{rl}GM& =2^{-59+60}\\ & \approx 2^1\\ & \approx 2\end{array}$$

似乎結果和第二個方法一樣，在沒有 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 的冪的運算，所以在浮點數編碼成二進位後，就只要做指數部分的運算了

$$\begin{array}{rl}GM& =[(M_0\cdot M_1\cdots \cdot M_n)\times (2^{{\alpha }_0}\cdot 2^{{\alpha }_1}\cdots 2^{{\alpha }_n}){]}^{\frac{1}{n}}\\ & =[(2^{{\beta }_0}\cdot 2^{{\beta }_1}\cdots 2^{{\beta }_n})\times (2^{{\alpha }_0}\cdot 2^{{\alpha }_1}\cdots 2^{{\alpha }_n}){]}^{\frac{1}{n}}\\ & =2^{\frac{1}{n}\sum _{i=0}^n({\alpha }_i+{\beta }_i)}\end{array}$$

浮點數的

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

$\beta$

$$\begin{array}{rl}& M=2^{\beta }\\ & \Rightarrow \beta =lo{g}_{2}M\end{array}$$

Exponent 的總結果

$E_{total}$ 會是

$$\begin{array}{rl}{E}_{total}& =\frac{1}{n}\sum _{i=0}^n{\alpha }_i+{\beta }_i\\ & =\frac{1}{n}\sum _{i=0}^n({\alpha }_i+lo{g}_2{M}_i)\\ & =\frac{1}{n}\sum _{i=0}^n[{\alpha }_i+lo{g}_2(frac\cdot 2^{-23}+1)]\\ & \approx \frac{1}{n}\sum _{i=0}^n{\alpha }_i\end{array}$$

接著會發現

$lo{g}_2(frac\cdot 2^{-23}+1)$ 的部分計算結果恆為零，因為

$frac$ 本來就只有 23 個位元，乘上

$2^{-23}$ 就相當於右移 23 個位元直接變成 0，所以只剩

$\frac{1}{n}\sum _{i=0}^n{\alpha }_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;
    
    e_sum = e_sum >= 0xff ? 0 : e_sum;

    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 中對浮點數進行二進位編碼，所以可以將

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

我們期望

$x$ 是

$k$ 倍的

$ln2$ 並且加上誤差

$r$，所以將 x 分解成

$$x=k\cdot ln2+r$$

所以

$k$ 應近似於

$$k=round(x\cdot \frac{1}{ln2})$$

$e^x$ 則可以分解成

$$e^x=e^{k\cdot ln2+r}=e^{k\cdot ln2}\times e^r=2^k\times e^r$$

$r$ 是

$x$ 和

$k\cdot ln2$ 的誤差，也就是一個正確值

$hi$ 與近似值

$lo$ 的差

$$r=x-k\cdot ln2=hi-lo$$

$$|r|<=0.5\times ln2\approx 0.34657359028$$

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

$e^r$ 部分

解

$e^r$ 和解

$e^x$ 的數值範圍並不一樣，這種拆解方式可以防止直接算

$e^x$ 時益位

以往計算

$e^x$ 都是直接使用泰勒展開式

$$\begin{array}{rl}{e}^r& =1+r+\frac{r^2}{2!}+\frac{r^3}{3!}+\frac{r^4}{4!}\cdots \\ & =1+r+\frac{r^2}{2}+\frac{r^3}{6}+\frac{r^4}{24}+\cdots \end{array}$$

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

$R(r^2)$

$$\begin{array}{rl}R(r^2)& =r\times \frac{(e^r+1)}{(e^r-1)}\\ & =r\times \frac{2+r+\frac{r^2}{2}+\frac{r^3}{6}+\frac{r^4}{24}+\cdots }{r+\frac{r^2}{2}+\frac{r^3}{6}+\frac{r^4}{24}+\cdots }\\ & =r\times \frac{\frac{2}{r}+1+\frac{r}{2}+\frac{r^2}{6}+\frac{r^3}{24}+\cdots }{1+\frac{r}{2}+\frac{r^2}{6}+\frac{r^3}{24}+\cdots }\end{array}$$

將

$u=\frac{r}{2}+\frac{r^2}{6}+\frac{r^3}{24}+\cdots$ 代入

$$R(r^2)=r\times \frac{\frac{2}{r}+1+u}{1+u}$$

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

• 第一項
  $\frac{2}{r}/(1+u)$

$$\frac{\frac{2}{r}}{1+u}=\frac{2}{r}\times (1-u+u^2-u^3+u^4+\cdots )$$

• 第二項
  $1/(1+u)$

$$\frac{1}{1+u}=1-u+u^2-u^3+u^4+\cdots$$

• 第三項
  $u/(1+u)$

$$\frac{u}{1+u}=u\times (1-u+u^2-u^3+u^4+\cdots )$$

所以

$R(r^2)$ 前 5 項是

$$\begin{array}{rl}R(r^2)& =r\times [(\frac{2}{r}+u+1)\cdot (1-u+u^2-u^3+u^4+\cdots )\\ & =2+\frac{r^2}{6}-\frac{r^4}{360}+\frac{r^6}{15120}-\frac{r^8}{604800}\cdots \end{array}$$

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

$2^{-59}$ 內

$$R(z)\approx 2+P_1\cdot z+P_2\cdot z^2+P_3\cdot z^3+P_4\cdot z^4+P_5\cdot z^5$$

,

$z=r^2$，而其中的 5 個係數就是

$$P1=\frac{1}{6}=1.66666666666666019037e-01$$

$$P2=-\frac{1}{360}=-6.61375632143793436117e-05$$

$$P3=\frac{1}{15120}=6.61375632143793436117e-05$$

$$P4=-\frac{1}{604800}=-1.65339022054652515390e-06$$

$$P5=4.13813679705723846039e-08$$

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

由

$R(r^2)$ 得到

$e^r$

$$\begin{array}{rl}& R(r^2)=r\times \frac{(e^r+1)}{(e^r-1)}\\ ⟹& e^r=\frac{r+R{(}^2)}{R(r^2)-r}\\ ⟹& e^r=1+\frac{2r}{R(r^2)-r}\end{array}$$

用

$R(z)$ 代入計算

$$e^r=1+\frac{2r}{R(z)-r}$$

此時討論如果

$r$ 接近 0:

$$R(r^2)\approx 2+\mathrm{極}\mathrm{小}\mathrm{值}$$

$$R(r^2)-r\approx 2-r+\mathrm{極}\mathrm{小}\mathrm{值}$$

計算

$\frac{2r}{2-r}$ 時，分母中較大的

$2$ 和較小的

$r$ 相減可能會失去精度，如:

假設支持 15 位的精準度，

$r=0.0001$:

$$2-0.0001+\frac{{0.0001}^2}{6}=1.99990001\bar{6}$$

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

$1.999900016666667$

所以整理成

$$e^r=1+r+\frac{r\cdot c(r)}{2-c(r)}$$

$$c(r)=r-(P_1\cdot r^2+P_2\cdot r^4+P_3\cdot r^6+P_4\cdot r^8+P_5\cdot r^{1}0)$$

因為

$c(r)$ 算起來會很小，所以

$2-c(r)$ 相對

$2-r$ 會更接近 2，減少精度損失，而且可以讓

$1+r$ 單獨計算

算出

$e^r$ 後，就可以用 my_ldexpf 算出

$e^x=2^k\times e^r$