浮點數運算基礎 (Floating-Point Fundamentals)
重點總覽
IEEE-754 二進位浮點標準自 1985 年被採用後，包含 CUDA 架構在內的主流運算系統都已實作此標準。浮點運算因精度有限，會做 rounding（捨入），使得「看似相同」的計算（不同結合順序、是否使用 FMA、是否沖掉 denormal）會得到不同位元結果，且常與 host CPU 結果不同——但這不代表 GPU 算錯。本篇涵蓋浮點格式、特殊值、FMA、CUDA 對 IEEE-754 與 C/C++ 的符合程度，以及內建算術運算子的精度與相關 nvcc 旗標。
| 項目 | 重點 |
|---|---|
| 為何 host/device 結果不同 | rounding、FMA contraction、加乘不滿足結合律、編譯器優化、library 實作差異 |
| 編碼三欄位 | sign（1 bit）、exponent（含 bias）、significand/mantissa |
| 預設 rounding | round-to-nearest-ties-to-even（rn）；CUDA 支援四種模式 |
| denormal（subnormal） | IEEE-754 必備，較慢；-ftz=true 可沖成 0（含於 --use_fast_math） |
| FMA | a*b+c 只 round 一次，更精確，可防 subtractive cancellation |
| IEEE-754 符合度 | 遵循 754-2019，但無動態 rounding mode、無例外偵測、SNaN 不 signaling |
+ - * / 精度 |
符合 IEEE-754，rn 下最大 0 ULP；-prec-div=false 後 /→2 ULP、1/x→1 ULP |
| 可移植預設旗標 | -ftz=false、-prec-div=true、-prec-sqrt=true，且不開 --use_fast_math |
Floating-Point Introduction
浮點格式與編碼
二進位浮點以三欄位編碼：
- Sign：1 bit，正負號。
- Exponent：以 base-2 指數，並偏移一個 bias。
- Significand（mantissa/fraction）：小數部分。
normal 值的數值公式：
(-1)^sign × 1.mantissa × 2^(exponent - bias)
subnormal（denormal）值改為：
(-1)^sign × 0.mantissa × 2^(1 - bias)
- 整數位
1.在 normal 值中是隱含的（implicit bit）。 - single bias = 127，double bias = 1023。
- 例：
-192 = (-1)^1 × 2^7 × 1.5，exponent 7 → float 編碼7+127=134，double →7+1023=1030。 - 因 fraction 位數有限，並非所有實數都能精確表示（如 2/3 = 0.10101010... 無限循環需 rounding）。
IEEE-754 定義的二進位格式對應 CUDA 型別：
| 位寬 | 別名 | CUDA/C++ 型別 |
|---|---|---|
| 16-bit | half-precision | __half |
| 32-bit | single-precision | float |
| 64-bit | double-precision | double |
| 128-bit | quad-precision | __float128 / _Float128 |
Normal 與 Subnormal
- exponent 欄位不全為 0 也不全為 1 → normal。
- 最小正 normal（
FLT_MIN）與 0 之間有一大段空隙，遠大於FLT_MIN到次小 normal 的間距。 - subnormal（denormal）：exponent 全 0、significand 至少一位非 0，用來填補此空隙，提供漸進的精度損失而非突然 round 到 0。
- denormal 是 IEEE-754 必備，但運算較昂貴。
不需嚴格精度時，可用 nvcc -ftz=true（flush-to-zero）關閉 denormal 提升效能；此旗標也包含在 --use_fast_math 內。
特殊值
| 特殊值 | 編碼 | 行為重點 |
|---|---|---|
| Zero | exponent 與 significand 全 0 | 有 +0 與 -0 兩種表示；+0 == -0 為 true |
| Infinity | exponent 全 1、significand 全 0（恰兩種編碼） | 飽和算術，溢位得 ±Infinity |
| NaN | exponent 全 1、significand 非全 0（共 2^(mantissa+1) − 2 種） | 任何含 NaN 的算術得 NaN |
- 不定型
Inf*0、Inf-Inf、Inf/Inf、0.0/0.0、sqrt(-1.0)皆得 NaN。 - 任何有序比較（
<、<=、>、>=、==）含 NaN 都回傳 false，包含NaN == NaN（非反身性）；只有NaN != NaN回傳 true。 - qNaN（quiet）傳播錯誤，significand 最高位=1；sNaN（signaling）用於引發 invalid 例外，最高位=0。確切 bit pattern 為 implementation-defined。可用
cuda::std::numeric_limits<T>::quiet_NaN/signaling_NaN取值。
Associativity（結合律不成立）
浮點加法與乘法不滿足結合律，因為每步可能需 rounding：
(A + B) + C 不一定等於 A + (B + C)
文件範例顯示：即使 A、B、C 全符合 IEEE-754，rn(rn(A+B)+C) 與 rn(A+rn(B+C)) 兩者結果互不相同，且其中 rn(A+rn(B+C)) 更接近數學真值。重點：運算順序影響 rounding error 的累積。
Fused Multiply-Add (FMA)
FMA 計算 a*b+c，只做一次 rounding（非 multiply 一次、add 一次共兩次），因此更精確。
- 防止 subtractive cancellation：乘法階段保留 double-width 乘積，即使加法時相消，乘積仍有足夠有效位。
- 範例（binary single）：
rn(A*A+B)得正確值，而分開算的rn(rn(A*A)+B)全部有效位被相消成 0。
FMA 對 NaN 傳播可能與兩步分開不同；多個 NaN 運算元時，payload 選擇順序為 implementation-defined。
CUDA 提供 FMA 的方式（float 與 double）：
| 方式 | 說明 |
|---|---|
x * y + z |
以 -fmad=true 或 --use_fast_math 編譯時 contract 成 FMA |
fma(x,y,z) / fmaf(x,y,z) |
C 標準函式庫 |
__fmaf_[rd,rn,ru,rz]、__fmaf_ieee_[...]、__fma_[rd,rn,ru,rz] |
CUDA 數學 intrinsic |
cuda::std::fma / cuda::std::fmaf |
CUDA C++ 標準函式庫 |
比對 CPU 與 GPU 結果時，須了解 host 平台是否啟用 FMA：GCC/Clang -mfma、NVC++ -Mfma、MSVC /fp:contract；硬體面如 x86 AVX2（-mavx2 / /arch:AVX2）、Arm64 Advanced SIMD (Neon)。host 是否 contract 會直接造成差異。
Dot Product 與 Rounding
同一個四元素內積，即使每步都符合 IEEE-754，不同實作策略仍可能得到略不同結果：
| 演算法 | 表示 |
|---|---|
| 1. 序列、乘加分開 | ((((a1·b1)+(a2·b2))+(a3·b3))+(a4·b4)) |
| 2. 序列、用 FMA | (a4·b4)+((a3·b3)+((a2·b2)+(a1·b1+0))) |
| 3. divide-and-conquer（parallel） | ((a1·b1)+(a2·b2)) + ((a3·b3)+(a4·b4)) |
IEEE-754 必備運算的跨實作保證：標準要求支援加、減、乘、除、sqrt、FMA、remainder、conversion、scaling、sign、comparison 等運算；對相同格式與相同 rounding mode，這些必備運算的結果保證在所有實作間一致。只有 IEEE-754 之外的 library 函式才不保證跨實作重現。
Rounding Modes：IEEE-754 定義四種，CUDA 全部支援，預設 rn；其他模式需用 intrinsic 對個別運算選擇。
| 模式 | 解讀 |
|---|---|
rn |
round to nearest, ties to even（預設） |
rz |
round towards zero |
ru |
round towards +∞ |
rd |
round towards −∞ |
影響 host/device 精度的因素
| 因素 | 重點 |
|---|---|
| Associativity | 加乘非結合，evaluation 順序影響 error 累積 |
| Fused Multiply-Add | 顯式呼叫或編譯旗標隱式啟用，依賴硬體支援 |
| Precision | 提高精度可減少 significance 損失，但 throughput 較低、用更多 register/記憶體 |
| Compiler Flags | -O3//O2 不改浮點語意，但 inlining、loop unrolling、向量化、CSE 可能改變結果；NVC++ 須 -Kieee -Mnofma 才完全 IEEE-754 |
| Library 實作 | IEEE-754 以外的函式不保證 correctly rounded，跨平台/架構可能不同 |
| Deterministic Results | 須相同硬體、編譯器版本/旗標、rounding mode/環境變數、輸入、thread 配置；atomic 排序依硬體排程，跑每次可能不同 |
善用 CUDA Math Libraries 與 C/C++ 標準數學函式：已最佳化且廣泛測試 edge case，避免手刻重造。
Floating-Point Data Types
CUDA 支援 Bfloat16、half、single、double、quad 精度，及供 tensor core 用的低精度格式（TF32、MX、4/6/8-bit）。
| 名稱 | 型別 | IEEE-754 | Header/Built-in | 需求 |
|---|---|---|---|---|
| Bfloat16 | __nv_bfloat16 |
✗ | <cuda_bf16.h> |
— |
| Half | __half |
✓ | <cuda_fp16.h> |
— |
| Single | float |
✓ | Built-in | — |
| Double | double |
✓ | Built-in | — |
| Quad | __float128/_Float128 |
✓ | Built-in；數學函式用 <crt/device_fp128_functions.h> |
host compiler 支援且 CC ≥ 10.0；拼法依 host compiler |
Bfloat16 需 Compute Capability 8.0 以上。TF32、microscaling (MX)、4/6/8-bit fp 等非通用計算用，主要服務 tensor core，詳見 CUDA Math API。
各型別範圍與精度（擇要保留關鍵數字）：
| 名稱 | 最大值 | 最小正 normal | 最小正 denormal | Epsilon |
|---|---|---|---|---|
| Bfloat16 | ≈ 3.39·10³⁸ | 2⁻¹²⁶ ≈ 1.18·10⁻³⁸ | 2⁻¹³³ | 2⁻⁷ |
| Half | 65504 | 2⁻¹⁴ ≈ 6.1·10⁻⁵ | 2⁻²⁴ | 2⁻¹⁰ |
| Single | ≈ 3.40·10³⁸ | 2⁻¹²⁶ ≈ 1.18·10⁻³⁸ | 2⁻¹⁴⁹ | 2⁻²³ |
| Double | ≈ 1.8·10³⁰⁸ | 2⁻¹⁰²² ≈ 2.22·10⁻³⁰⁸ | 2⁻¹⁰⁷⁴ | 2⁻⁵² |
| Quad | ≈ 1.19·10⁴⁹³² | 2⁻¹⁶³⁸² ≈ 3.36·10⁻⁴⁹³² | 2⁻¹⁶⁴⁹⁴ | 2⁻¹¹² |
- Bfloat16 與 Single 範圍幾乎相同（同 8-bit exponent），但 Bfloat16 mantissa 較少（Epsilon 2⁻⁷ vs 2⁻²³），精度低、動態範圍大。
- 用
cuda::std::numeric_limits（<cuda/std/limits>）查詢各型別屬性與範圍，含 MX 格式。 - 複數：C++ 標準庫
cuda::std::complex（<cuda/std/complex>）；CUDA 另提供cuComplex/cuDoubleComplex（cuComplex.h）。
CUDA and IEEE-754 Compliance
GPU 遵循 IEEE 754-2019，但有以下限制：
- 無動態可設 rounding mode；多數運算支援多個常數 IEEE rounding mode，透過特定命名的 device intrinsic 選擇。
- 無浮點例外偵測機制：所有運算如同 IEEE 例外永遠 masked。SNaN 編碼雖支援，但不 signaling，當作 quiet 例外處理。
- 浮點運算可能改變輸入 NaN 的 payload bit；abs、negation 等也可能不符 IEEE-754，NaN 的 sign 可能 implementation-defined 地被更新。
最大化結果可移植性，建議用 nvcc 預設旗標 -ftz=false、-prec-div=true、-prec-sqrt=true，且不開 --use_fast_math。注意：浮點 re-association 與 contraction（FMA）預設允許（如同 --fmad=true）。
浮點轉整數且 round 後落在目標整數範圍外時，IEEE-754 與 C/C++ 標準未明定。GPU 的 clamping 行為見 PTX ISA conversion 指令；但若非直接經 PTX 指令觸發，編譯器可能援引未指定行為 → undefined behavior、無效程式（如 __double2int_rz()）。可能與 host 不同。
Atomic 函式的 denormal 行為（不論 -ftz 設定）：
| Atomic 運算 | denormal 行為 |
|---|---|
| global memory 上的 atomic single-precision add | 永遠 flush-to-zero（等同 PTX add.rn.ftz.f32） |
| shared memory 上的 atomic single-precision add | 永遠支援 denormal（等同 PTX add.rn.f32） |
CUDA and C/C++ Compliance
浮點例外：與 host 不同，device 端數學運算子/函式不設 errno、不回報浮點例外。需要錯誤診斷時，使用者須自行對輸入輸出做 screening。
Undefined Behavior（UB）：
- 一般 UB：未初始化的浮點變數、超出生命週期、signed integer overflow、解參考無效指標。
- 浮點特有 UB：把浮點轉成無法表示的整數型別（含 NaN、infinity）→ UB。
- 浮點除以 0 不是 UB，而是 implementation-specific：符合 IEC-60559 (IEEE-754) 的實作（含 CUDA）產生 infinity，且不受編譯器優化影響。
- 無效浮點運算（如
0/0、Inf/Inf）產生 NaN，這是合法行為，不可誤判成 UB。 - 整數除以 0 才是 UB。
浮點字面值可移植性：C99/C++17 支援十六進位浮點字面值（base-2 精確），十進位字面值可能無法以 base-2 表示。兩者皆 round 到最近可表示值（implementation-defined），host 與 device 可能不同。
float f1 = 0.5f; // 十進位字面值，0.5 可精確表示
float f2 = 0x1p-1f; // 十六進位字面值，亦為 0.5
float f3 = 0.1f; // 0.1 無法以 base-2 精確表示
- run-time 評估受選定 rounding mode、FMA contraction、reassociation 設定、浮點例外影響（CUDA 不支援例外、預設 rn）。
- 編譯器對常數表達式可能用更高精度內部表示，並做 constant folding/propagation、CSE，導致最終值或比較結果不同。
C 標準數學庫：host 實作以平台相關方式對映 host libm；host 沒有的函式在 crt/math_functions.h 實作（如 erfinv()）；少見函式（如 rhypot()、cyl_bessel_i0()）僅 device 可用。host 與 device 實作獨立。
Floating-Point Functionality Exposure
CUDA 透過多種方式暴露數學功能，各有支援型別、host/device 可用性、是否受浮點旗標影響：
| 暴露方式 | 主要型別 | Host | Device | 受 FP 旗標影響（僅 float/double） |
|---|---|---|---|---|
| Built-in C/C++ 算術運算子 | float, double, __half, __nv_bfloat16, __float128/_Float128, cuda::std::complex |
✓ | ✓ | ✓ |
CUDA C++ 標準庫（<cuda/std/cmath>, cuda::std::） |
同上 + fp8/fp6/fp4 查詢* | ✓ | ✓ | ✓ |
CUDA C 標準庫 (Math API)（<math.h> 子集，免額外 header） |
float, double（half/bf16/fp128 有限） | ✓ | ✓ | ✓ |
| Non-standard CUDA 數學函式 (Math API) | float, double（half/bf16/fp128 有限） | 視函式而定（per-function） | ✓ | ✓ |
| Intrinsic 數學函式 (Math API) | float, double | ✗ | ✓ | 僅 -ftz=true |
* CUDA C++ 標準庫對 __nv_fp8_e4m3/e5m2/e8m0、__nv_fp6_e2m3/e3m2、__nv_fp4_e2m1 等小型別支援查詢（numeric_limits<T>、fpclassify()、isfinite()、isnormal()、isinf()、isnan()）。
Intrinsic 函式：較快但較不精確，只在 device 可用。其行為不受 -prec-div=false、-prec-sqrt=false、-fmad=true 影響；唯一例外是 -ftz=true（含於 --use_fast_math）。
__half/__nv_bfloat16需<cuda_fp16.h>/<cuda_bf16.h>。__float128/_Float128依 host compiler 與 CC；相關 Math API 函式需<crt/device_fp128_functions.h>且僅 device 可用。- 函式精度可能 host 與 device 不同（標準庫常底層呼叫 CUDA Math API）。
- 精度量化用 ULP（Unit in the Last Place）。
Built-In Arithmetic Operators
內建運算子 x+y、x-y、x*y、x/y、x++、x--、倒數 1/x（single/double/quad）符合 IEEE-754：rn 模式下保證最大 0 ULP 誤差，host 與 device 皆可用。
-fmad=true（含於 --use_fast_math）啟用乘加 contraction，對 float 的最大 ULP 影響：
| 旗標/運算 | 結果 | 最大 ULP（float） |
|---|---|---|
-fmad=true：x*y+z |
→ __fmaf_rn(x,y,z) |
0 ULP |
-prec-div=false：x/y |
→ __fdividef(x,y) |
2 ULP |
-prec-div=false：1/x |
— | 1 ULP |
預設 -prec-div=true 時 / 是 correctly rounded（0 ULP）；開 --use_fast_math（內含 -prec-div=false）後 / 與 1/x 變快但引入 1~2 ULP 誤差。
考試/測驗重點
| 題型 | 關鍵答案 |
|---|---|
| GPU 結果與 CPU 不同代表 GPU 錯了？ | 否；差異源自 rounding/FMA/結合順序，需謹慎解讀 |
| 浮點預設 rounding mode？ | round-to-nearest-ties-to-even（rn） |
| CUDA 支援幾種 rounding mode、怎麼選非預設？ | 四種（rn/rz/ru/rd）；用 intrinsic 對個別運算選 |
| GPU 能動態切換 rounding mode 嗎？ | 不能；只能透過特定命名的 device intrinsic 用常數模式 |
| denormal 是什麼、如何關閉？ | exponent 全 0 的 subnormal；-ftz=true（flush-to-zero，含於 --use_fast_math） |
| FMA 為何更精確？ | a*b+c 只 round 一次，並用 double-width 乘積防 subtractive cancellation |
啟用 x*y+z contract 成 FMA 的旗標？ |
-fmad=true（含於 --use_fast_math） |
| 浮點加乘滿足結合律嗎？ | 不滿足；順序影響 rounding error 累積 |
+0 == -0？NaN == NaN？ |
+0==-0 為 true；NaN==NaN 為 false（NaN!=NaN 為 true） |
| 哪些運算得 NaN？ | 0/0、Inf-Inf、Inf/Inf、Inf*0、sqrt(-1.0) |
| 浮點除以 0 是 UB 嗎？ | 否，implementation-specific，IEEE-754 實作得 infinity；整數除 0 才是 UB |
| CUDA 符合哪版 IEEE-754、有何限制？ | 754-2019；無動態 rounding、無例外偵測、SNaN 不 signaling |
| 哪些運算保證跨實作可重現？ | IEEE-754 必備運算（+、−、×、÷、sqrt、FMA、remainder、conversion、scaling、sign、comparison）在相同格式與 rounding mode 下保證一致；IEEE-754 之外的 library 函式不保證 |
| 最佳可移植性的 nvcc 旗標？ | -ftz=false、-prec-div=true、-prec-sqrt=true，不開 --use_fast_math |
| 採用可移植預設旗標後 FMA 會被關閉嗎？ | 不會；re-association 與 FMA contraction 預設仍允許（等同 --fmad=true），要關需另設 -fmad=false / NVC++ -Mnofma |
| global vs shared 上 atomic float add 的 denormal？ | global 永遠 flush-to-zero；shared 永遠支援 denormal（不論 -ftz） |
+ - * / 1/x 在 rn 下 ULP 誤差？ |
0 ULP（IEEE-754 correctly rounded） |
-prec-div=false 後 / 與 1/x 的 ULP？ |
x/y→__fdividef 2 ULP；1/x 1 ULP |
| Intrinsic 函式受哪些 FP 旗標影響？ | 不受 -prec-div/-prec-sqrt/-fmad 影響；唯 -ftz=true 例外 |
| device 端數學函式會設 errno 嗎？ | 否；不設 errno、不回報浮點例外 |
| Bfloat16 vs Single 的範圍與精度？ | 範圍幾乎相同（同 8-bit exponent）；Bfloat16 mantissa 少、精度低（Eps 2⁻⁷） |
| Bfloat16 需要哪個 CC？ | Compute Capability 8.0 以上 |
| 哪些型別是 IEEE-754？ | half、float、double、quad 是；Bfloat16 不是 |
| 精度量化單位？ | ULP（Unit in the Last Place） |