FMA 只捨入一次——先把整精度乘積算完，再捨進最終結果，比先乘後加省下一次誤差。作者照著一篇有 Coq 形式證明的論文寫了一份軟體實作，測試全部綠燈，直到另一位貢獻者拿真實硬體的答案一比對，冒出一組連他自己都對不上的輸入。
寫一個 FMA，順手撿到 stdlib 的 bug
Fused multiply-add——FMA——算的是 a×b+c，只捨入一次，不是先乘一次捨一次、再加一次捨一次：「Fused multiply-add (FMA) computes a * b + c with only one rounding error instead of two.」硬體裡通常直接有一條指令，但不是所有機器都有——Firefox 的硬體遙測顯示，「15% of machines in the Firefox hardware survey don't have AVX2 and hardware FMA that comes with it」，換算下來每七台裝置就有一台得靠軟體算。在這之前，Rust 的 std::simd 本來就有一份後援：「On machines without hardware FMA, Rust's std::simd gives up and runs scalar FMA on each f32 in [f32; 4] individually, which is slow.」——遇到沒有硬體 FMA 的機器，直接放棄向量化，一個一個 f32 慢慢算。shnatsel 想做的是更好的版本：「I wanted to do better in fearless_simd and provide an actually vectorized implementation.」底子是 Sylvie Boldo 與 Guillaume Melquiond 2008 年那篇《Emulation of FMA and correctly-rounded sums》，演算法本身在 Coq 裡有形式證明，理論上捨入邏輯不該出錯。先乘後加、各自捨入一次，跟一次算完再捨入一次，兩者在大多數輸入上答案一樣；但當乘積剛好卡在捨入邊界附近，中間多的那一次捨入，有機會把答案推到跟正確答案差一格的地方——這就是為什麼 FMA 要求只捨入一次，而不是把兩個現成的浮點運算疊在一起用。這類底層數值函式一旦進了 libm 或 std，就會被無數上層程式碼直接信任——不會有人每次呼叫 fmaf() 之前自己再驗一次答案對不對。
我的浮點乘加，一路綠燈到硬體比對那一刻
在浮點函式庫的測試慣例裡，「已知值」通常指的是那些特別容易踩雷的邊界——正負零、無窮大、NaN、剛好卡在指數進位那一格的數字；shnatsel 沒有在文章裡列出他挑的每一個已知值，只說「some known values plus some random tests」。FMA 作用在三個 f32 上，可能的輸入組合有多少種，shnatsel 自己算過：「the total number of possible inputs is 2 to the 96th power. It would take the world's largest supercomputer only 500 years to try them all」——窮舉測試在時間尺度上直接出局。剩下能做的，只有挑幾個已知邊界值，加上大量隨機測試：「random tests that try a million subnormals, and random tests that try a million values that should require fixup」，subnormal 一百萬次、需要捨入修正的情況再一百萬次——浮點捨入實作裡常見的做法，是用 guard、round、sticky 這幾個額外位元記住「被捨掉的那段裡還有沒有東西」，交給最後一步決定怎麼捨入；這一步在 normal 數值上有現成的位元技巧可以抄近路，換到 subnormal 身上，因為少了那個隱含的 1，同樣的技巧就會失靈。CI 也沒閒著，fearless_simd 的建置設定會在每一種 SIMD 等級的模擬器裡再跑一輪：「CI configured to run tests in an emulator for every supported SIMD level, on top of running them on the host」。這一整套測試全部通過。
點任一卡片看這層測試擋不擋得住這顆 bug · 4 層
四層測試，哪一層真的擋住了這顆 bug？
detail
L1 · 形式證明的演算法基礎
證明保證的是紙上那個演算法本身——合理的推測是，它不保證某一份被抄進 C 或 Rust 的具體位元判斷式，每一步都跟論文完全對應。
L2 · 已知值＋百萬級隨機測試
兩百萬次隨機測試全部通過，卻沒抽到那一組會踩雷的輸入——bug 藏得比隨機取樣的覆蓋率窄。
L3 · CI 跨 SIMD 等級模擬
合理的推測是，這一層測的是同一份程式碼在不同 SIMD 路徑之間彼此一致，而不是拿去跟另一份獨立實作或硬體結果做外部比對。
L4 · 跟真實硬體結果差分比對
前三層都沒攔下的東西，這一步一次就撞上——真正抓到 bug 的是外部差分比對，不是任何一層自己寫的測試。
通過歸通過，跟「正確」是兩回事。真正讓答案現形的不是這幾層測試，是另一位開發者 @awxkee 在 fearless_simd 的 PR 底下留言，貼出一組跟真實硬體 FMA 算出來的結果對不上的輸入——shnatsel 自己的說法是「@awxkee appeared and posted some inputs on which my implementation diverged from the hardware results」。百萬級隨機測試、跨 SIMD 等級的模擬 CI，都沒撞到這一組輸入；一次跟硬體的差分比對就撞上了。這種「拿另一個獨立來源做差分比對」的做法，在浮點數正確性驗證上特別關鍵——浮點捨入的正確性沒有辦法只靠讀程式碼確認邏輯來驗證，得真的把可能的路徑跑過一次，跟一個信得過的參照值比對。已知值測試涵蓋的是想得到的邊界；隨機測試涵蓋的是統計意義上的覆蓋率；跟硬體或另一份獨立實作做差分測試，涵蓋的是想不到、但確實存在的角落——三種方法各自守住不同的縫，缺一種，就會在那個縫上失守。
假設一：是不是我抄的證明演算法本身有洞？
形式證明能做的事，其實跟單元測試不是同一個量級——它證的是「對所有可能輸入，這個演算法都會給出正確答案」，不是「我測過的這些輸入都對」。這正是為什麼 shnatsel 一開始選擇照著論文寫，而不是自己從頭設計捨入邏輯：省下的不只是時間，是把「邏輯本身是否正確」這件事的驗證責任，轉嫁給已經做過的人。問題出在轉嫁之後那一步——把證明過的演算法，一個字元一個字元地翻譯成能跑的程式碼，這段翻譯過程本身沒有形式保證。第一個念頭是演算法本身：Boldo 與 Melquiond 那篇論文是形式證明過的，如果連 Coq 都點頭，捨入邏輯照理該是對的。但形式證明保證的是紙上那個 correctly roundedFMA 的定義本身——先當作在無限精度下算完 a×b+c，只捨入這最後一步，不是先乘一次捨一次、再加一次捨一次。 演算法，合理的推測是，它沒有辦法保證某一份被抄進 C 或 Rust 的具體位元判斷式跟論文的每一步都對得上。真正出問題的，是演算法轉譯成程式碼之後，處理 halfway case捨入到最接近偶數時，剛好卡在兩個可表示浮點數正中間的輸入——原始判斷式的注解「/* not a halfway case */」講的就是這種情況。 的一行判斷式，而這行判斷式只認得 subnormal指數欄位全部是 0、沒有隱含前導 1 的浮點數——數值小到已經沒辦法再用正常的「指數＋尾數」結構去逼近。 以外的數值。換句話說，問題不是出在被證明過的那一步，而是出在證明範圍之外、負責判斷現在是不是卡在正中間的那一小段位元運算——這一小段本身沒有被形式驗證覆蓋到。
假設二：會不會只是隨機測試的樣本數不夠？
如果懷疑是測試樣本不夠多，得先弄清楚 bug 藏得有多窄。原始的反例是三個位元組：a = 0x97000800、b = 0x1cfff001、c = 0x00010002，正確答案應該是 0x00010001，軟體 FMA 卻算出 0x00010002——兩個答案在整數位元上只差 1，也就是差 1 個 ULP（unit in the last place，浮點數在目前這個位置能表示的最小間距，中間沒有任何一個浮點數卡得進去）。shnatsel 給的誤差量級是「off by about 15 parts per million of the correct result」。把 c 的 hex 展開來看，指數欄位剛好是全 0，也就是說這組反例裡的加數 c 本身就落在 subnormal 範圍——這是從已知的位元組直接算出來的，原文沒有另外點出這件事。把兩個答案的尾數展開看，正確答案的尾數是 65537，bug 答案的尾數是 65538，正好差 1——相對誤差是 1÷65537，換算成百萬分之一大約是 15，跟 shnatsel 給的「about 15 parts per million」對得上。這種需要特殊處理的捨入情況，本來就不常見：「you hit it less than once in a million when processing values in the [-1, 1) range」。百萬次隨機 subnormal、百萬次隨機 fixup，命中率終究是隨機的；而這顆 bug 只需要落在 1 個 ULP 寬的縫裡就會現形。把這個縫放回原本的空間看更清楚：2 的 96 次方種輸入裡，只有極小一段會踩到這道 1 個 ULP 寬的界線——百萬次等級的隨機抽樣，命中率趨近於零，不是因為測試寫得不夠好，而是統計覆蓋率天生贏不了這種等級的稀疏度。
拖動滑桿，看正確答案跟 bug 答案的位元距離 · 17 個相鄰值
真正的答案：一行只認得 normal 數值的判斷式
真正的問題出在一行判斷式：if ((u.i & 0x1fffffff) != 0x10000000) /* not a halfway case */。compiler-builtins 是 Rust 標準函式庫背後、提供底層數值運算（像是浮點四則運算、整數除法）的那個 crate——fmaf 這類函式最終就是從這裡編譯進 std。它的 issue #1262 裡，shnatsel 把話講得很白：「This detection of points exactly halfway between representable ones only works for normal values, but not for subnormals.」IEEE 754 的 float32 用 1 個符號位元、8 個指數位元、23 個尾數位元表示一個數。normal 數值的尾數前面有一個隱含的 1——不占位元、但永遠存在，公式是 (-1)^sign × 1.mantissa × 2^(exponent-127)。指數欄位全部歸零的時候，這個隱含的 1 消失，數值改用 (-1)^sign × 0.mantissa × 2^(-126) 表示——這就是 subnormal，用來讓浮點數在逼近零的時候平滑地失去精度，而不是直接斷崖式跳到零。少了那個隱含的 1，任何假設「尾數前面一定有個 1」的位元技巧，套到 subnormal 頭上都會算錯——這正是那行判斷式踩到的坑：它只認得 normal 數值的位元排列，subnormal 數值指數欄位固定是 0，內部結構不一樣，同一個遮罩套上去，判斷結果就不對。
這顆 bug 不是 Rust 一家的事。Rust 的 std::simd／compiler-builtins fmaf 是這樣，musl libc 的 fmaf() 也是同一顆——musl 的郵件論壇在 8 月 10 日收到一模一樣的位元反例（x = 0x97000800、y = 0x1cfff001、z = 0x00010002，musl 算出 0x00010002，硬體 FMA 是 0x00010001）。musl 那邊的回報者具名 Sergey Davidoff——從 bug 內容跟時間點來看，合理推測就是 shnatsel 本人。兩邊的實作往上追，都指向同一份出處：musl 的原始碼把這段 fmaf() 的實作來源記成 FreeBSD。
| 實作／專案 | 平台範圍 | 回報 | 目前狀態 |
|---|---|---|---|
| Rust std::simd／compiler-builtins fmaf | 缺硬體 FMA 的 x86、32 位元 ARM 等平台 | issue #1262（8/9）、PR #1270 | 修正已送出，still awaiting review |
| musl libc fmaf() | 同上，musl 常見於嵌入式、Alpine 等環境 | mailing list 回報（8/10） | minimal fix 與 substantial rewrite 都已提出，尚未合併 |
| FreeBSD（原始出處） | musl 沿用至今 | 由 musl 原始碼回溯認出出處 | 原文沒有另外提到 FreeBSD 自身是否處理 |
三份實作共用同一行判斷式。shnatsel 指出 musl 那份檔案把實作歸給 FreeBSD 時，標的版權年份是「with copyright years 2005-2011」——版權年份標的是主張權利的年份，不等於程式碼寫成的時間，但這段邏輯被 musl 沿用到今天是可以確定的。原因不難理解：c = 0x00010002 這種輸入要同時滿足兩件事才會踩雷——加數本身要落在 subnormal 範圍，又剛好卡在捨入的 halfway 邊界上。b 是 normal 數值，指數落在 −70 附近；c 卻是 subnormal，指數欄位全 0，實際大小落在 2⁻¹³³ 這個量級——一個運算元在「正常」範圍、一個貼著浮點數能表示的地板，這種組合在真實工作負載裡本來就少見。三個運算元裡，a 本身也不是普通數字——它是負數，指數同樣落在很小的量級，但仍然是 normal；整組反例湊在一起，剛好覆蓋了「正常範圍的小數字乘上正常範圍的小數字，加上貼著地板的 subnormal」這種特定組合。這也是這類 bug 特別頑固的原因——同一段捨入邏輯，靠著程式碼複製，從一份實作散布到另一份，中間沒有人重新驗證過。真正拿去跟硬體 FMA 比對過的，目前只有 shnatsel 這份 fearless_simd 實作；合理的推測是，musl 跟它上游的 FreeBSD 版本，在這之前都沒有經過同一等級的差分測試，才會讓同一顆 bug 安穩地活了十幾年。對還在用 32 位元 ARM 或缺 AVX2 的 x86 平台跑數值運算的人來說，這不是一顆需要立刻停機搶修的 bug——百萬分之十五等級的相對誤差，大多數應用感覺不到。但如果程式碼剛好依賴同一段浮點運算在不同機器上得出逐位元一致的結果，像是某些物理模擬、重播測試、或是拿雜湊值比對計算結果，這顆 bug 就有機會讓人在排查「為什麼這兩台機器結果不一樣」時繞不少彎路。
拖動手把，比較 b、c 在整個浮點數量級光譜上的距離 · 範圍 2⁻¹⁵⁰ 到 2⁻⁶⁰
水平軸是浮點數的量級（2 的次方，非線性壓縮）。左邊灰色段是 subnormal 範圍，指數欄位固定是 0；跨過 2⁻¹²⁶ 這條分界才進入 normal 範圍。c 落在灰色段裡，b 落在灰色段右邊很遠的地方——兩者相差超過 60 個 2 的次方。
這顆 bug 會咬到誰，又還沒補上多久
波及的硬體範圍不算小，但也不是全部：「On x86 this is only a problem on cheap and/or very old Intel without AVX2, and on Chinese Hygon x86 chips that don't have FMA at all. 64-bit ARM is unaffected.」64 位元 ARM 沒事——原文只給了這個結論，沒有交代原因；真正的重災區是 32 位元 ARM 這類舊架構——shnatsel 自己的判斷是「this will probably haunt 32-bit ARM and various embedded systems for years to come」，但這句話他自己加了「probably」，也承認沒有把每一份嵌入式工具鏈的組合語言都讀過一遍去確認：「I sure am not going to read that much ARM assembly」，他另外補了一句：「I would not be surprised if the buggy FreeBSD implementation got copy-pasted into lots of different toolchains for various obscure platforms.」
修 bug 的過程還意外撞見第二顆、完全無關的 bug。Rust 維護者 tgross35 在測試 f128 這個實驗型別時，發現 32 位元 ARM hardfloat 平台上有 calling convention 對不上的問題：「On arm32 we use aapcs_on_arm for all float intrinsics so f16/f32/f64 get passed in GPRs. However, it looks like on hardfloat targets LLVM is expecting them to be passed in float registers.」同一段程式碼路徑，兩顆完全不相干的 bug——一顆潛伏十幾年的 subnormal 捨入 bug，一顆全新型別上路才浮現的 ABI 不匹配，都不是靠讀程式碼讀出來的，而是靠實際跑起來、拿去跟另一個獨立來源比對才現形。這顆 ABI 不匹配已經透過 PR #1292 關閉；shnatsel 自己送出的 FMA 修正 PR #1270，撰稿當下的狀態是「It's still awaiting review.」他在準備這份 PR 的過程裡，另外跑了一次 LLM 檢查自己有沒有引入新的迴歸——「It occurred to me to run an LLM to search for bugs I could have introduced, and it found a counter-example」，真的抓到一個會弄壞 floating-point exception flags 的反例。他也補了一句：「using LLMs for analysis is permitted by the Rust LLM policy, under certain conditions」，算是先打了預防針——在一篇靠人工排查揪出 stdlib bug 的文章裡，特別註明這句話本身也是一種交代方法的態度。
拖動手把回顧回報順序 · 6 個事件
這顆 bug 安靜到沒有任何錯誤訊息——fmaf() 不會 panic、不會回傳 NaN，只是把答案往旁邊挪一格。這種看起來完全正常、只是差一點的失效模式，比崩潰難抓得多：崩潰會逼人去查，一個差 15 ppm 的答案通常直接被當成正確結果用掉。至於這一切在真實世界裡有多大殺傷力，shnatsel 自己給的答案很直接：「I have no idea.」他能想到的具體場景是「deterministic simulations running across different machines will diverge」——同一段模擬程式碼在不同機器上跑，只要有一台落到軟體 FMA 路徑、又剛好踩中這個誤差範圍，結果就對不齊。musl 那邊，「a minimal fix and a substantial rewrite」都已經提出、也經過其他開發者的意見來回修過，但「it's still not merged as I'm writing this」。一顆潛伏在 FreeBSD 程式碼裡十幾年沒被抓到的 bug，離「補完」還有一段距離。「提出修法」跟「修法真的合併進 release」之間那段距離，在開源基礎建設裡並不罕見，尤其是這種底層、影響範圍難以估計、又需要多方 review 才敢動的程式碼。對還在用受影響平台的人來說，與其等上游合併，不如先照 issue #1262 跟 musl mailing list 給的位元反例，自己驗一次目前用的版本。
Next time：形式證明保證的是紙上的演算法，不保證那份被抄進 C 或 Rust 的位元判斷式跟著證明走；百萬次隨機測試贏過的是統計覆蓋率，不是跟硬體或另一份獨立實作的差分比對——真正抓到這顆 bug 的，是一次跟硬體結果的比對，不是任何一層自己寫的測試。