Linux 核心專題: 浮點數運算案例探討

執行人: dcciou
解說影片1
解說影片2
GitHub

Reviewed by Daniel-0224

A common choice for fixed-point arithmetic in Linux kernel development is a scaling factor of 1024 (corresponds to a shift of 10 bits), which means that
there will be 22 bits for the integer part and 10 bits for the fractional part.
This setup provides a resolution of 2^-10 = 0.00097656, which sufficiently meets the precision requirements of many kernel-related tasks without the need for
extremely fine scales

可以標示一下此段文章的出處。

已修正，感謝指正。

Reviewed by chloe0919

1. hackmd 撰寫請符合中文和英文字元之間要有空白字元

透過右移操作縮減val，右移的位數是local_n 除以LOAD_AVG_PERIOD 的結果。

可否說明為何右移的位數是 local_n 除以LOAD_AVG_PERIOD 的結果

已修正，感謝指正

Reviewed by fennecJ

• repo 中提供的測試 程式碼 以及本筆記中摘錄的程式碼含有中文註解，請改為使用美式英語撰寫。

• 可否在你的筆記中針對 kahan summation, pairwise summation 的計算方法進行簡介，並嘗試分析這兩種計算法有何特性以及各自適合的運算場景？

• 你的執行程式 runtimetest_fixed_point 在 user space 下進行計時，可能導致結果受到 interrupt/preemption 影響而有所誤差。可以透過撰寫 kernel module 關閉 interrupt/preemption 取得更精準的結果。

已修正，感謝指正及指教

Reviewed by popo8712

在處理浮點數乘 10 的實作中，您提到嘗試 Kahan summation 和 pairwise summation 來消除誤差，但未能成功。能否詳細解釋為什麼這些方法沒有達到預期效果？是否考慮過其他消除浮點數運算誤差的方法？

根據本人對浮點數的了解，這是因為單精度浮點數的格式限制了有效位數，超過小數點後六位仍然無法避免誤差，要同時滿足如此精度以及維持可表示的數值範圍應該是無解；但若是數值不大，或許可以使用定點數，以小數點後六位的精度來說，需要平移至少 20 bits。

Reviewed by mesohandsome

筆記中的圖片全都無法看到，需要修改權限。

圖片權限應該無法修改，可能請您重整看看。

Reviewed by dockyu

我們可以把* 10變成* 8 * 1.25，也就是

$2^3\ast (1+1/2)$，在程式碼操作上使用加﹑除﹑位元操作來完成，程式碼如下

0.25 應該是

$\frac{1}{4}$ ，應該改為

$2^3\ast (1+\frac{1}{4})$

已修正，感謝指正

Reviewed by marsh-fish

建議註解用英文書寫。

感謝提醒

Reviewed by lintin528

在你的 normal -> overflow 或其他的情況中，可以嘗試自己設計一些補償機制並訂定資料集來觀察。

了解，感謝指教

Reviewed by david965154

• 直接乘10:

Time used:0.000025 secs

• bitwise 8+2:

Time used:0.000101

• fixed point:

Time used:0.000021

那關於誤差的分析呢，單看這這樣的數據時間似乎差的並不多，可以提供一下在你的測試資料下精度為何

精度部分直接乘與bitwise皆為小數點後六位，定點數部分為小數點後三位。

任務簡介

回顧 IEEE 754，探究其原理並進行案例探討，接著研究 Linux 核心原始程式碼的定點數運算，理解其設計和實作考量。

TODO: IEEE 754 背景介紹

務必依循 IEEE 754 規格書來探討浮點數標準和關鍵考量

浮點數的類型

• 根據 IEEE754(2008 ver) ch3.1 提到，浮點數的format共分為五種：Three binary formats, with encodings in lengths of 32, 64, and 128 bits, and two decimal formats, with encodings in lengths of 64 and 128 bits。
• 另外，在 ch3.2 提到，浮點數共有三種值的類型：Triples, infinity, NaN

Triples

Triples指的是：sign, exponent, significand。

• sign : 決定數值的正負
• exponent : 決定數值的規模
• significand : 提供數值的精確度
  具體呈現方式如下:
  

修正圖片的存取權限。

當 e=emin and 0<m<1，則歸類為 Subnormal numbers

Infinity(6.1)

對無限操作數的運算通常是精確的，因此不會觸發異常(exception)，除了某些例外

• addition(∞, x), addition(x, ∞), subtraction(∞, x), or subtraction(x, ∞), for finite x
• multiplication(∞, x) or multiplication(x, ∞) for finite or infinite x ≠ 0
• division(∞, x) or division(x, ∞) for finite x
• squareRoot(+∞)
• remainder(x, ∞) for finite normal x
• conversion of an infinity into the same infinity in another format.

例外

• ∞ 為無效的操作數
• ∞ 從有效操作數的 overflow 產生 (see 7.4) 或 /0 (see 7.3)
• remainder(subnormal, ∞) signals underflow

使用 LaTeX 符號

Not a number(6.2)

NaN有兩種： signaling and quiet

• signaling NaNs: 用於表示一些特殊情況，例如未初始化的變數或一些超出 IEEE 754 標準範圍的操作（如複數中的無窮大或極寬的數值範圍）
• quiet NaNs: 表示未定義或無法計算的結果，如:0/0﹑無限大。由於 quiet NaNs 不會跳警告而是繼續執行，因此在每次 quiet NaNs 出現的時候做紀錄將會有助於 debug，這部分 IEEE 讓 user 決定如何紀錄診斷資訊。為了便於傳播 NaN 中包含的診斷資訊，應盡可能多地將該資訊保留在 NaN 操作結果中。除了 fusedMultiplyAdd 操作外，任何涉及到 quiet NaNs 的操作皆不會觸發異常。

浮點數的規格

在 IEEE 754 ch3.6 可以看到浮點數的規定如下:

將常用的 binary32, binary64 圖像化表示如下:

數據範圍如下:

修正圖片存取權限

浮點數的計算

根據 IEEE 754(2008 ver) ch4.3　提到，在計算時一律先計算至無限精度及範圍，再根據條件做四捨五入

Except where stated otherwise, every operation shall be performed as if it first produced an intermediate
result correct to infinite precision and with unbounded range, and then rounded that result according to one
of the attributes in this clause.

捨入的方法

由於浮點數並不像實數具有連續性，也不像自然數般均勻分布，而是像電子能階一樣的不均勻分布，因此數值的捨入方式變得需要因地制宜。這邊一共有兩個類型，共五種方案

Attributes to nearest

這邊討論的是尚未碰到最大值的處理方案，一個浮點數的最大值可以表示成：

$b^{emax}(b-({\displaystyle \frac{1}{2}}{b}^{1-p}))$

• roundTiesToEven: 當結果不準確 (precise) 且有兩個同等近的浮點數可供選擇時，選擇最後一位為偶數的那個。二進制的預設捨入法。
• roundTiesToAway: 當結果不準確且有兩個同等近的浮點數可供選擇時，選擇具有較大絕對值的，也就是負更負，正更正。

Directed rounding attributes

這邊討論的是用戶可以選擇的捨入方向

• roundTowardPositive: 當進行捨入時，結果總是會向正無窮方向捨入，也就是 無條件進位。應用場合可以是計算獲利的時候，避免客戶有所損失。
• roundTowardNegative: 向負無窮捨入是將結果捨入為比實際計算結果小的最接近的可表示數值，也就是 無條件捨去。應用場合可以是預算編列的時候，避免因為捨入問題超過上限。
• roundTowardZero: 結果被捨入到零的方向，正數結果向下捨入，負數結果向上捨入，使得結果的絕對值總是小於或等於實際計算結果的絕對值。也就是正數無條件捨去，負數無條件進位。應用場合會是在想要避免變號的情況。

注意：

The default rounding-direction attribute for results in decimal formats is language defined, but should be roundTiesToEven.

結果

根據 IEEE754(2008 ver) ch5.1 ，符合本規定的浮點數操作產生的結果大致可分成四類

一般計算操作(General-computational operations)

產生浮點或整數結果，並對所有結果進行四捨五入，有可能觸發浮點異常。

安靜計算操作(Quiet-computational operations)

產生浮點結果，且不會觸發浮點異常訊號。本操作又可再分成兩類(5.5)。

• Sign bit operations: 它們只影響正負位元。這些操作對浮點數和 NaN 一視同仁，且不觸發訊號異常。這些操作可能會傳播非規範 (non-canonical) 編碼。
  • sourceFormat copy(source): 複製操作，將浮點數x 複製到相同格式的另一個位置，不改變正負位元。
  • sourceFormat negate(source): 取反操作，將浮點數 x 的正負位元反轉，等同於執行 x * -1，但更有效率，因為只改變正負位元。
  • sourceFormat abs(source): 絕對值操作，將浮點數 x 的正負位元元設為正（0），即取得 x 的絕對值。
  • sourceFormat copySign(source, source): 複製符號操作，將浮點數 x 複製到相同格式的另一個位置，但使用 y 的正負位元。
• Decimal re-encoding operations: 與十進位格式相關的轉換操作，使得浮點資料可以在十進位和二進位編碼之間轉換。
  • decimalEncoding encodeDecimal(decimal): 將十進位浮點數編碼成為特定的十進位編碼形式。
  • decimal decodeDecimal(decimalEncoding): 將十進位編碼解碼回十進位浮點數。
  • binaryEncoding encodeBinary(decimal): 將十進位浮點數編碼為二進制浮點數的形式。
  • decimal decodeBinary(binaryEncoding): 將二進位編碼的浮點數解碼為十進位浮點數。

訊號計算操作(Signaling-computational operations)

不會產生浮點結果，可能會觸發異常;比較操作屬於此類(5.6)。

非計算操作(Non-computational operations)

不會生成浮點結果，也不會觸發異常(5.7)。

此外，根據結果格式和操作數格式之間的關係，操作也可以分成兩類

• homogeneous operations: 浮點操作數和浮點結果都具有相同的格式
• formatOf operations: 指示結果的格式，與操作數的格式無關

異常(exception)

在 IEEE754(2008ver) ch7.1 提到，異常共有五種。這些事件通常是由於運算超出了正常的數值範圍或邏輯，需要特殊處理。

註: flag 不同於 signal ， flag 是一種狀態指示，不會讓程式運行中止。

Invalid operation　(7.2)

當沒有有用的可定義結果時，才會發出無效操作異常信號。對於以浮點格式生成結果的操作，發出無效操作異常的預設結果應為　quiet NaN，該 NaN 應提供一些診斷資訊（見 6.2）。這些操作包括：

• 任何包含 signal NaN 的通用計算或訊號計算，一些特定的轉換除外。
• 乘法運算：0 乘無窮大（∞）或無窮大乘0。
• 融合乘加運算（fusedMultiplyAdd）：0 乘以∞ 再加上一個常數 c，或∞ 乘以0 再加上 c，除非 c 是 quiet NaN。 如果 c 是 quiet NaN，則是否觸發異常取決於具體實作。
• 加法或減法或融合乘加：無窮大的減法，例如+∞ 加上−∞。
• 除法操作：0 除以0 或無限大除以無限大。
• 餘數運算：當 y 為零或 x 為無窮大且兩者皆非NaN 時所進行的餘數運算。
• 開平方根運算：運算元小於零時嘗試開平方根。
• 量化運算（quantize）：結果不符合目標格式，或一個運算元是有限的而另一個是無限大。

對於未以浮點格式生成結果的操作，發出無效操作異常信號的操作包括：

• 浮點數轉換成整數：當來源浮點數是 NaN、無限大（infinity），或轉換結果超出整數格式可表示的範圍時，轉換操作會觸發無效操作異常。
• 當操作數未排序時，通過表 5.2 中列出的無序信令謂詞進行比較
• 當對數函數 logB 的輸入是 NaN、無限大或零，且結果的預期格式是整數時，操作會觸發無效操作異常。

Division by zero (7.3)

當一個運算在有限操作數上執行，並且該運算的結果被定義為無限大時。 這通常涉及到除數為零而被除數為非零的有限數。預設結果如下:

• 對於除法：當除數為零且被除數是一個有限的非零數時，無限大的正負是兩個運算元的 XOR 結果。也就是：如果被除數和除數的符號相同（都是正或都是負），結果為正無窮（+∞）；如果符號不同（一個正一個負），則結果為負無窮（−∞）。
• 對於logB(0)：當logB 函數的輸入為0 且輸出格式為浮點格式時，結果的符號為負（−∞）。 原因是對數函數在數學上對零取對數是未定義的，而按照常見的數學規定，log(0) 被認為趨向於負無窮。

Overflow (7.4)

當在沒有指數範圍限制的情況下，經過捨去的浮點運算結果的絕對值超過了目標格式能表示的最大有限數值。以單精度來說，若指數部分的二進位表示法超過 11111111，則觸發 overflow 。捨入方式如第４章所提:

• roundTiesToEven﹑roundTiesToAway：所有溢出的結果都被捨入到無窮大(∞)，符號與中間結果（即捨去前的結果）相同。
• roundTowardZero：所有溢出的結果都被捨入到該格式能表示的最大有限數，符號與中間結果相同。
• roundTowardNegative：正溢位被捨入到該格式的最大有限數，負溢位被捨去到負無窮（−∞）。
• roundTowardPositive：負溢位被捨入到該格式的最負的有限數，正溢位被捨去到正無窮（+∞）。

Underflow (7.5)

當檢測到微小的非零結果時，應觸發 underflow exception 。對於二進位格式，應在下列兩種方式中選擇一種做為檢測方法：

• 捨入後的 underflow ：當非零結果在沒有指數範圍限制下，經舍入後的結果介於
  $\pm b^{emin}$ 之間（bemin 是可表示的最小非零正規數的指數，範圍不包含端點，即嚴格介於的概念）。
• 捨入前的 underflow ：當非零結果，如果沒有指數範圍和精確度限制，計算的結果介於
  $\pm b^{emin}$ 之間。

預設的異常處理方式如下：

• 無論檢測微小性的方法如何，總是會交出一個舍入後的結果，可能是零、次正規數，或
  $\pm b^{emin}$。
• 如果舍入結果不精確：即舍入的結果和沒有指數範圍和精度限制時的計算結果不同，那麼將會觸發 underflow flag ，並觸發 inexact exception。
• 如果舍入結果精確：舍入結果與沒有指數範圍和精度限制時的計算結果相同，則不會觸發任何flag，也不會 觸發 　inexact exception。

Inexact

除非另有說明，否則當運算的四捨五入結果不精確時（即，它與指數範圍和精度均無界時計算的結果不同），應發出不精確異常的訊號。四捨五入或溢出的結果應交付到目的地。

TODO: 針對單精度浮點數運算，提出利用 bitwise 運算達到「乘 10」的操作

不能用 mul

一般情況

也就是 i/o 都是 normal numbers 的情況。
我們可以把* 10變成* 8 * 1.25，也就是

$2^3\ast (1+1/4)$，在程式碼操作上使用加﹑除﹑位元操作來完成，程式碼如下

注意書寫規範：

• 使用 lab0 規範的程式碼書寫風格，務必用 clang-format 確認一致
• 注意看程式碼規範的「clang-format 工具和一致的程式撰寫風格」，筆記列出的程式碼依循指定的風格。

#include <stdio.h>
#include <stdint.h>

float multiplyByTen(float num) {
    // Step 1: x + x/4
    float result = num + num / 4.0;

    // Step 2: Multiply by 8 using bitwise manipulation
    uint32_t* bits = (uint32_t*) &result;
    uint32_t exponent = (*bits >> 23) & 0xFF;  // take exponent part

    if (exponent == 0 || exponent == 0xFF) {
        // handle special condition：zero，Inf， NaN
        return result;
    }

    // bit shift：*8 = shift 3 bits
    int newExponent = (int)exponent + 3;

    if (newExponent >= 0xFF) {
        // overflow can be set as Inf
        *bits = (*bits & 0x80000000) | 0x7F800000;
    } else {
        // update exponent part
        *bits = (*bits & 0x807FFFFF) | (newExponent << 23);
    }

    return result;
}

int main() {
    float value = 1.234567f;
    float result = multiplyByTen(value);
    printf("Original: %f\n", value);
    printf("Multiplied by 10: %f\n", result);
    return 0;
}

從浮點數的規範可以推知，單精度浮點數的有效位數約為七位，因此測試值的輸出大機率會出現誤差。輸出結果如下。

Original: 1.234567
Multiplied by 10: 12.345671

測試少一位是否有誤(正確):

Enter a floating-point number: 1.23456
Original: 1.234560
Multiplied by 10: 12.345600

很明顯的可以看到誤差出現，透過 IEEE-754 Floating Point Converter可以看到實際被存入的數值與輸入不同。

修改圖片的存取權限

再利用另一個 converter得到的結果顯示，該輸入為 inexact ，因此我開始嘗試消除誤差。

• 嘗試 1: kahan ， 沒用

Enter a floating-point number: 1.234567
Original: 1.234567
Multiplied by 10 using Kahan sum: 12.345671

• 嘗試 2: Pairwise Summation ，失敗

Enter a floating-point number: 1.234567
Original: 1.234567
Multiplied by 10 using pairwise summation: 12.345671

• 嘗試 3: 計算過程採 double ，儲存時取 single ，失敗

Enter a floating-point number: 1.234567
Original: 1.234567
Multiplied by 10 using double precision: 12.345671

• 解決辦法 1:改 double。輸出正確，但會消耗更多儲存空間。

Enter a floating-point number: 1.234567
Original: 1.234567
Multiplied by 10: 12.345670

• Q: 0* 10
  A: 很明顯的 0* 10被歸類在 normal number

Enter a floating-point number: 0
Original: 0.000000
Multiplied by 10: 000000

normal -> overflow

(signaling NaN)
bit shift 無解，就是 overflow
Q: 有沒有辦法以同樣位元數來表示結果
A: 可以，但尾數的間距會變大導致誤差，需要額外補償機制

改進你的漢語表達。

overflow -> overflow

bit shift 無解，繼續 overflow

Original: inf
Multiplied by 10: inf

underflow -> normal

應該有解，會出現quiet nan

underflow -> underflow

或許可以bit shift 到 normal，之後再想要怎麼還原(多乘的倍數如何校正)

TODO: 研究 Linux 核心原始程式碼的定點數運算

參照《Demystifying the Linux CPU Scheduler》關於定點數運算的描述，以 loadavg 和 PELT 等 Linux 核心原始程式碼裡頭的實作，探究其如何處理定點數的除法。
分析能否善用上述利用 bitwise 運算達到「乘 10」的操作」達到簡化運算成本。

定點數概述

將數值縮放後儲存，根據定義的小數點位置決定縮放的比例

常見格式選擇

在老師所著的 Demystifying the Linux CPU Scheduler 之第69頁提到，小數部分通常放大

$2^{10}$ 倍，原因是符合 Linux 所需的精度(約為小數點後三位)。整數的部份用 22 bits 儲存，小數為 10 bits。

A common choice for fixed-point arithmetic in Linux kernel development is a scaling factor of 1024 (corresponds to a shift of 10 bits), which means that
there will be 22 bits for the integer part and 10 bits for the fractional part.
This setup provides a resolution of 2^-10 = 0.00097656, which sufficiently meets the precision requirements of many kernel-related tasks without the need for
extremely fine scales

範例如圖: 圖中儲存的數值為 123.25，其二進位為1111011.01，位移10 bits 後存進去

操作方式

轉換

int f = (int) (v * (1 << fixedpoint_shift));

乘法

int f = (f1 * f2) >> fixedpoint_shift;

除法

int f = (f1 / f2) = (f1 * (1 / f2)) >> fixedpoint_shift;

PELT(Per-entity load tracking)

概述

為改善進程的CPU時間公平性以及提升效率而被提出的方法，首次出現於 Linux3.8 中。
本法按照每個調度實體（可以是單一進程或整個控制組）來追蹤負載，其中每個實體的負載是基於它們的可運行時間來計算，包括實際運作的時間和等待運作的時間。系統將這些時間劃分為1024 微秒的周期，追蹤每個週期中實體的狀態變化。

• 實體的負載貢獻：在每個週期中，如果實體是可運作的（無論它是否真正佔用了CPU），都被視為對系統負載有貢獻。 這種方式不僅考慮了實體實際佔用 CPU 的時間，還包括了它在就緒佇列中等待的時間。

• 負載衰減演算法：為確保負載追蹤的資料不會因時間推移而變得不相關，核心引入了衰減因子（y），用於減少歷史負載資料的影響。 具體來說，32 毫秒前的負載貢獻權重為目前的一半，越早的負載貢獻對目前的調度決策影響越小，這個概念類似一階馬可夫鏈。
  

  

定點數部分

static u64 decay_load(u64 val, u64 n)
{
    unsigned int local_n;

    if (unlikely(n > LOAD_AVG_PERIOD * 63))
        return 0;

    local_n = n;

    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); 
    return val;
}

透過右移操作縮減val，右移的位數是local_n 除以LOAD_AVG_PERIOD 的結果。

Loadavg

概述

衡量系統負載情況的一個重要指標，它表示系統在特定時間間隔內的平均活躍進程數。 活躍進程指的是那些處於運作狀態或等待狀態（即正在使用或等待使用CPU 資源）的進程。

定點數部分

以除法為例，這邊是透過bit shift來完成的，這部分的程式碼如下:

static unsigned long
fixed_power_int(unsigned long x, unsigned int frac_bits, unsigned int n)
{
    unsigned long result = 1UL << frac_bits;  // set initial number

    if (n) {
        for (;;) {
            if (n & 1) {
                result *= x;
                result += 1UL << (frac_bits - 1);  // rounding
                result >>= frac_bits;              // devided by bit shift
            }
            n >>= 1;
            if (!n)
                break;
            x *= x;
            x += 1UL << (frac_bits - 1);  // rounding
            x >>= frac_bits;              // down accruacy loss by shit bits
        }
    }

    return result;
}

成本比較

以下皆為測試 10000 次的運行時間，可以看到定點數運算比直接乘來的快20%左右，但精度只能到小數點後三位。 bitwise可能是我程式碼不夠簡潔才導致運算時間大幅提升。

• 直接乘 10:

Time used:0.000025 secs

• bitwise 8+2:

Time used:0.000101

• fixed point:

Time used:0.000021

詳述實驗方法和解讀。