2025-02-{18,20} 課堂問答簡記

編輯本頁，紀錄課堂問答內容並留下你的 GitHub 帳號

jserv

連鋼琴師都用 Linux

comment: 人家是 NiceChord (好和弦)，我們是 NiceCode (好程式)

devarajabc

課程簡介和注意須知第 21 頁提及的程式碼:

#include <stdio.h>
int main() {
    float sum = 0.0f;
    for (int i = 0; i < 10000; i++) sum += i + 1;
    printf("Sum: %f\n", sum);
    return 0;
}

其執行結果是 Sum: 50002896.000000，與 1 + 2 + 3 + … + 10000 實數加總的預期結果 50005000 不同，該如何解釋？以 IEEE 754 標準，針對浮點數表示和對應的位元逐一說明。

為什麼會這樣呢？
因為並不是每個整數都能以單精度浮點數的型態表示！
舉例來說，將單精度浮點數 A 做加法並輸出：

float A = 16777216.0;
printf("%f \n",A);
printf("%f \n",A+1);
printf("%f \n",A+2);
printf("%f \n",A+3);
printf("%f \n",A+4);

會得到錯誤的結果，16777217 和 16777219 消失了

16777216.000000 
16777216.000000 
16777218.000000 
16777220.000000 
16777220.000000

Why ? 讓我們來觀察浮點數轉換的過程：

image

整數 \(16777217\) 的二進位制表示式為 \(1000000000000000000000001_2\)(25 位元)
\(16777217 = 1.000000000000000000000001_2* 2^{24}\)
其中 \(000000000000000000000001_2\)(24 位元) 會被存入 mantissa field, 不過因為單精度浮點數的 mantissa field 只有 23 位元，所以實際存入的只有 \(00000000000000000000000_2\)(23 位元)
如此一來，\(16777217\) 和 \(16777216\) 以單精度浮點數型態來表示都會是 \(1.0* 2^{24}\)

難道超過 \(16777216\) 的奇數以單精度浮點數來表示就會自動減一嗎？

錯！大錯特錯，觀察下面的例子：
整數 \(16777219\) 的二進位表示式為 \(1000000000000000000000011_2\) (25 位元) ，最後一位元的 1 被捨棄了，以單精度浮點數來表示就是 \((1.0+2^{-23}) \times 2^{24}=16777218\) 怎麼會這樣呢？ 為什麼實際輸出的結果是 \(16777220\) 呢？

(我不知道，四捨五入的規則是什麼呢？Rounding to even ?)
感謝 Dennis40816 的說明
在規格書 4.3.1 中提到 roundTiesToEven：

The floating-point number nearest to the infinitely precise result shall be delivered; if the two nearest floating-point numbers bracketing an unrepresentable infinitely precise result are equally near, the one with an even least significant digit shall be delivered.

無法以浮點數型態表示的數字會以最接近其數值的數來表示，若「最接近的數字」有兩個，則會以 LSB 為 0 的數值來表示。注意，此處的 "least significant digit" 指的是 mantissa field 內的 least significant digit.

那問題來了，如何知道「最接近的數字」是多少呢？

16777216 -> 16777216
16777217 -> 16777216 (-1)
16777218 -> 16777218 
16777219 -> 16777220 (+1)
16777220 -> 16777220
16777221 -> 16777220 (-1)
16777222 -> 16777222
16777223 -> 16777224 (+1)

根據輸出的結果，我們可以觀察到：

  • 整數範圍 [\(2^{24}\), \(2^{25}-1\)] 中只有偶數能以浮點數型態表示
  • mantissa field 的 LSB 為 0 的數會出現三次
  • 以浮點數型態表示的數值和實際數值的差異存在 [0, -1, 0, +1] 的循環

接下來我們來討論 26 位元的整數以浮點數型態表示的形況:
\(33554433\) 以二進位表示為 \(10000000000000000000000001_2\)
其 significand 有 25 位，末端兩位元無法表示。換句話說，只有四的倍數能以浮點數型態表示(末兩位為零) ，配合 roundTiesToEven ，可以得到下表：

33554432 -> 33554432 
33554433 -> 33554432 (-1)
33554434 -> 33554432 (-2)
33554435 -> 33554436 (+1)
33554436 -> 33554436 
33554437 -> 33554436 (-1)
33554438 -> 33554440 (+2)
33554439 -> 33554440 (+1)
33554440 -> 33554440 
33554441 -> 33554440 (-1)
33554442 -> 33554440 (-2)
33554443 -> 33554444 (+1)
33554444 -> 33554444 
33554445 -> 33554444 (-1)
33554446 -> 33554448 (+2)
33554447 -> 33554448 (+1)
33554448 -> 33554448 
33554449 -> 33554448 (-1)
33554450 -> 33554448 (-2)
33554451 -> 33554452 (+1)
33554452 -> 33554452 
33554453 -> 33554452 (-1)
33554454 -> 33554456 (+2)
33554455 -> 33554456 (+1)
33554456 -> 33554456

根據輸出的結果，我們可以觀察到：

  • 整數範圍 [\(2^{25}\), \(2^{26}-1\)] 中只有四的倍數可以以浮點數形式表示
  • 誤差存在循環的情況，即 [0, -1, -2, +1, 0, -1, +2, +1]

既然知道「誤差值」會循環出現，是否可以利用餘數來計算誤差值呢？

int diff24[4] = {0, -1, 0, 1};
int diff25[8] = {0, -1, -2, 1, 0, -1, 2, 1};

以 sum%4 作為 diff24 的索引值
但考量到 sum 本身可能就是被錯誤近似的值
採用 sum%4 + i%4 作為索引值
問題： 如果 sum 無法表示呢？

wx200010

課程簡介和注意須知第 22 頁提及的程式碼:

int main() {
    float sum = 0.0f, corr = 0.0f; /* corrective value for rounding error */
    for (int i = 0; i < 10000; i++) {
        float y = (i + 1) - corr; /* add the correction to specific item */
        float t = sum + y; /* bits might be lost */
        corr = (t - sum) - y; /* recover lost bits */
        sum = t;
    }
    printf("Sum: %f\n", sum);
    return 0;
}

Q1: 英語的 "correction" 和 "corrective" 是什麼意思？(查閱英英字典，得知其用法)

[correction(名詞)](https://dictionary.cambridge.org/zht/%E8%A9%9E%E5%85%B8/%E8%8B%B1%E8%AA%9E-%E6%BC%A2%E8%AA%9E-%E7%B9%81%E9%AB%94/correction)：改正；修改；更正。 [corrective(形容詞)](https://dictionary.cambridge.org/zht/%E8%A9%9E%E5%85%B8/%E8%8B%B1%E8%AA%9E-%E6%BC%A2%E8%AA%9E-%E7%B9%81%E9%AB%94/corrective)：改正的，糾正的。

漢語在描述概念的分野時，很難有精準的詞彙，因此你務必查閱 English-English dictionary。例如 instruction 和 command 是英語沒有任何共通詞源的單字，但因為漢語分別翻譯為「指令」及「命令」，結果導致一堆人根本無法區分二者，本例可看出現代漢語的侷限。

由於先前直接查閱中文翻譯實在很難看出這兩者的差異，我重新查閱了其他免費的英英字典：

  • correction(noun)
    1. a change that makes something more accurate than it was before
    2. the act or process of correcting something
  • corrective(adjective)
    1. designed to make something right that was wrong before

對於 correction 的英譯，著重於「使某事物比之前更準確的一項更改」或「一個能修正某事物的行為或過程」；而 corrective 的英譯，則是「旨在糾正先前錯誤的某事物」，能用於形容「某個用來糾正錯誤的東西」。

為了更好地體會這兩者在語意上的不同，我們來回頭看看程式碼註解吧！
在程式碼中，有兩行註解使用到這兩項單字，分別是：

float sum = 0.0f, corr = 0.0f; /* corrective value for rounding error */

float y = (i + 1) - corr; /* add the correction to specific item */

整份程式碼在執行過程中會修正浮點數運算因精度不足造成的 rounding error，變數corr 的作用便是計算出缺失值，並在後續的計算中將缺失值加回結果，來讓結果更準確。

在第一段程式碼旁，原註解是 /* corrective value for rounding error */，意思是「用於處理rounding error的修正值」，主要是對變數corr的存在意義進行說明。
若將原註解改為/* correction for rounding error */，意思就會變成「修正 rounding error」，然而這行程式碼只是在宣告變數，並未實際修正rounding error，因此語意上會不太恰當。

而在第二段程式碼旁，原註解是/* add the correction to specific item */，意思是「對特定項目進行修正」，強調了這一行程式碼的主要目的。
若將原註解改成/* add the corrective value to specific item */，意思會變成「將修正值加入特定項目」，雖然的確是這行程式碼所執行的操作，但這樣會變成強調「加入某個數值」的操作，而不是「對某個東西進行修正」的主要目的，所以重點會變不一樣。

Q2: Kahan summation algorithm 的原理？限制是什麼？

由於討論前需要先了解 rounding error 的知識，這邊先簡單說明 rounding error 為何會發生。
首先，單精度浮點數分成 sign (1 bit), exponent (8 bits), fraction (23 bits) 來表示，fraction (23 bits) 用於表示二進位的小數部分，而在小數部分前面會有個隱藏的 1 (即 1.xxxxxxxxxx...)。
當時我透過這些概念做了份猜想：「一個單精度浮點數能夠儲存的精度可能是由fraction的23 bits+前面隱藏的1 bit掌管」。下面的程式範例則可以讓大家感受看看這份猜想：

#include <stdio.h>
int main()
{
    int nums[4] = {0x400001, 0x800001, 0x1000001, 0x8000007};
    float f;
    for(int i = 0; i < 4; ++i){
        f = nums[i];
        printf("nums[%d] = %ld, f = %.2f\n", i, nums[i], f);
    }
}
/* output : 
nums[0] = 4194305, f = 4194305.00
nums[1] = 8388609, f = 8388609.00
nums[2] = 16777217, f = 16777216.00     
nums[3] = 134217735, f = 134217728.00   
*/

我們發現，整數0x1000001與0x8000007在轉換成單精度浮點數會有誤差產生。
若將這兩個整數改以二進位來觀察，會發現最高位為1的bit位置與最低位1的bit位置間隔皆超過了23個bits，因此若將這些整數轉換成單精度浮點數(float)，會因精度限制(24 bits)導致只能留住從最高位為1開始的24個bits，剩下的那些低位bits則會被丟棄。

image

註：這張示意圖只是我目前對浮點數轉換上的理解，不確定是不是完全正確的。

關於捨去低位 bits 的部分，這篇文章 提到
Decimal places: The exact or closest value of the decimal places of the storage factor.
另外我自己用 FloatConverter 玩也是一樣的結果（會看後面那段要被捨棄的 bits 第一位是 0 或 1 ，如果 0 就捨去 1 就進位，感覺就是把十進位的四捨五入改成二進位零捨一入xd），也因此 rounding error 變大邊小都有可能，以 1 加到 100000 為例，得出的結果會是 5000050176.000000 ，較正確值 5000050000 還大 moonjam

回頭來看課程簡介和注意須知第 22 頁的部分程式碼:

float sum = 0.0f, corr = 0.0f; /* corrective value for rounding error */
for (int i = 0; i < 10000; i++) {
    float y = (i + 1) - corr; /* add the correction to specific item */
    float t = sum + y; /* bits might be lost */
    corr = (t - sum) - y; /* recover lost bits */
    sum = t;
}

y = (i + 1) - corr 是這輪要加進結果的值 (第一次發生 rounding error 前 corr 皆是 0)，而t = sum + y;則是這輪運算後所得到的實際結果值。通常來說，t應與我們理想的結果值 sum + y 相同，然而在發生 rounding error 時，t就會與原先預期的結果值sum + y不同。因此透過corr = (t - sum) - y的運算，便能將這輪因 rounding error 導致的誤差計算出來，從而能夠校正下次的計算。

為什麼 corr = (t - sum) - y 能夠計算出這輪的誤差值呢？
y 是這輪預期加進結果的值，而 (t-sum) 則是這輪實際加進結果的值，將兩者相減後，corr便得到這輪因 rounding error 而造成的誤差值了，便能在下次的計算中去校正結果值。

要注意的是，若改成 corr = t - (sum + y)，會因為(sum + y)本身會觸發 rounding error 的關係，導致兩者相減後會恆為0而變得毫無意義，因此將 t - (sum + y)改為 (t - sum) - y 也是為了避免計算誤差值時又再次發生 rounding error 的窘境。

就很像是(L + R) / 2 可以改寫為L + (R - L) / 2，來避免(L + R)發生overflow的問題
數學上等價的算式，對有限的二進位來說可能會有不同的結果。

目前的校正方式什麼時候會失敗呢？
那便是使y = (i + 1) - corr發生rounding error的情況，y的缺失值就會直接消失，造成不準確。
例如當i不低於0x1000000時，(i + 1)的部分就可能會開始發生rounding error了，很危險！
在使用這種校正方法前，必須要考慮清楚或明確限制輸入值的範圍，否則會無濟於事。

對Kahan summation algorithm的閱讀還未完成，因此內容可能會再更新。

\(\to\) Mado: 完全使用定點數避開浮點數運算的進階視窗系統

LinMarc1210

Q: 如果在 for 迴圈的最後一輪發生 corr 不為 0，那少加的值會不會沒被加進 sum

會，實際測試從 1 加到 9997，使用投影片上的程式碼會算出 49975004.000000 而非正確答案 49975003，然而即便在迴圈外再做 sum -= corr，輸出的答案仍然會是 49975004.000000，因為 49975003 無法用 float 表示，應該只能分別輸出 sum 和 corr 或將兩者轉為 int 再相加才能解決吧？ moonjam

devarajabc

Q: Kahan summation algorithm 要如何處理那些「無法被浮點數型態表示」的數字?，例如 \(50005001\)
修改投影片第 22 頁的程式：

-    float sum = 0.0f, corr = 0.0f; /* corrective value for rounding error */
+    float sum = 1.0f, corr = 0.0f; /* corrective value for rounding error */

預期輸出應為 \(50005001\)，但實際輸出為 \(50005000\)

moonjam

Q: 正版的 IEEE 754 規格書 好像要錢，有什麼免費且不侵害智慧財產取得的方式嗎？

用 "Draft" 來找，例如 IEEE Standard for Floating-Point
Arithmetic

Select a repo