IEEE 754 Floating Point

浮點數標準

大名鼎鼎的 IEEE 754 是定義 32-bit 浮點數的官方文件，於 1985 年由發布第一個版本。
下圖^[1] 顯示 Float32 各位元的功能:

Image Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

位元可以劃分為三部分:

sign bit (藍色): 31^th 位，1 表示負，共 1 位
exponent bits (綠色) : 23^th - 30^th 位，共 8 位
fraction bits (紅色): 0^th - 22^th 位，共 23 位

計算的方法如下:

(- 1)^{s} \times (1 + \frac{f r a c t i o n}{2^{23}}) \times 2^{(e - 127)}

浮點數具有幾種特殊值:

exp 為 0b1111_1111 && fraction 為 0 && sign 為 1: -INF
exp 為 0b1111_1111 && fraction 為 0 && sign 為 0: INF
exp 為 0b1111_1111 && fraction 為非 0 && sign 為 0: NaN

要將 10 進位小數轉換成浮點數格式需要兩個步驟

10 進位 → 2 進位
2 進位 → 浮點數表達 (sign, exponent, fraction)

10 進位數值是如何轉換成 2 進制呢?

假設有

X

為一 10 進制數字，可分為整數

N

和小數

F

兩部分:

X = N + F

N

使用除 2 取餘可得，餘數順序分別是

b_{0}, b_{1}, \dots

，當被除數小於 2 時停止

N = b_{k} \cdot 2^{k} + b_{k - 1} \cdot 2^{k - 1} + \dots + b_{1} \cdot 2^{1} + b_{0} \cdot 2^{0}

F

使用乘 2 取整方式 (乘 2 後看整數是否為 1) 得到，當被乘數等於 0 時停止

F = b_{- 1} \cdot 2^{- 1} + b_{- 2} \cdot 2^{- 2} + \dots

讓我們以 9.73 為例

\begin{aligned} N & = 9 = 2^{3} + 2^{0} = 1001_{b} \\ F & = {0.73}_{10} \\ \approx {0.101110101110000101001}_{b} \\ X & = 1001_{b} + {0.101110101110000101001}_{b} \\ = 1.001_101110101110000101001_{b} \cdot 2^{3} \end{aligned}

將

X

換成以浮點數表示:

sign:
$0$
exp:
$127 + 3 = 130 = 2^{7} + 2^{1} = 1000_0010_{b}$
fraction:
$0011_0111_0101_1100_0010_100 (1)_{b}$ ，最後一位被省略

組合起來就是:

\begin{matrix} \underset{sign}{\underset{⏟}{0}} \underset{exponent}{\underset{⏟}{1000_0010}} \underset{fraction}{\underset{⏟}{0011_0111_0101_1100_0010_100}} \end{matrix}

透過印出 9.73 float 的 binary format 得到同樣的結果:

float f = 9.73;
unsigned int* pu32_f = (unsigned int*)&f;
unsigned int u32_f = *pu32_f;
printf("[%d] ", (u32_f) & (1 << 31), 1, 0);
for (int bit = 30; bit >= 0; bit--)
{
    int c = 30 - bit;
    printf("%d", (u32_f) & (1 << bit) ? 1 : 0);
    if ((c - 3)% 4 == 0)
    {
      printf("_");
    }
}

浮點數的表達極限在哪裡?

由於不可能透過有限位數表達無限多的實數，因此透過位元表達勢必存在極限，下圖^[3]很好地說明符點數的分佈會隨著數字增加，而使最小精度單位 (Unit in Last Place, ULP) 以 2 的指數倍下降。例如:

[2^{23}, 2^{24})

時 ULP 是 1，而對於

[2^{24}, 2^{25})

則是 2。

Image Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

浮點數特別的極限:

當
$X = 2^{23}$ : 由於整數部份已經佔據 23 bits，因此無任何小數資訊能被保存，ULP 為 1，若繼續增加則 ULP 會變成 2, 4, …

範例觀察
接下來讓我們用程式觀察幾個 ULP 計算範例:

import numpy as np

num = np.float32(123456789.0)  # 轉換為 float32
next_up = np.nextafter(num, float('inf'), dtype=np.float32)  # 找比它大的最近浮點數
next_down = np.nextafter(num, float('-inf'), dtype=np.float32)  # 找比它小的最近浮點數

print("="*40)
print(f"  num      = {num}")
print(f"  next_up  = {next_up}")
print(f"  next_down= {next_down}")
print(f"  ULP 誤差範圍 = {next_up - num}")
print("="*40)

# 比較不同的數值，顯示結果更清楚
test_values = [
    (123456795.99, 123456800),
    (123456796.0,  123456800),
    (123456803.99, 123456800),
    (123456804.0,  123456800),
    (123456804.1,  123456800),
    (123456804.0,  123456796.0),
]

print("\n浮點數比較結果:")
print("-"*40)
for val1, val2 in test_values:
    res = np.float32(val1) == np.float32(val2)
    print(f"  np.float32({val1}) == np.float32({val2}) → {res}")
print("-"*40)

========================================
  num      = 123456792.0
  next_up  = 123456800.0
  next_down= 123456784.0
  ULP 誤差範圍 = 8.0
========================================

浮點數比較結果:
----------------------------------------
  np.float32(123456795.99) == np.float32(123456800) → False
  np.float32(123456796.0) == np.float32(123456800) → True
  np.float32(123456803.99) == np.float32(123456800) → True
  np.float32(123456804.0) == np.float32(123456800) → True
  np.float32(123456804.1) == np.float32(123456800) → False
  np.float32(123456804.0) == np.float32(123456796.0) → True
----------------------------------------

觀察以下兩組比較:

針對兩鄰近的浮點數表示 123456792.0 和 123456784.0，由於 123456789.0 更靠近前者，因此會選用前者的表示方式
針對兩鄰近的浮點數表示 123456800.0 和 123456808.0，由於 123456804.0 與兩者同樣接近，最後選擇了前者，Why?

這跟浮點數的捨入方法有關:

當兩邊距離相同時，捨入至最低有效位 (LSB) 為偶數方 (IEEE754-2019, 4.3.1 roundTiesToEven)

123456800.0 的 LSB 是 0 而，12345808.0 則是 1，故捨入前者

123456800.0 的 IEEE 754 表示: 01001100111010110111100110100100 (0x4CEB79A4)
123456808.0 的 IEEE 754 表示: 01001100111010110111100110100101 (0x4CEB79A5)
123456800.0 的 LSB: 0
123456808.0 的 LSB: 1

能正確被浮點數表示的數字稱機器可表示數 (machine-representable number, MRN)，以 MRN 為圓心，半徑為 0.5*ULP 範圍內的數字都會以 MRN 紀錄，剛好在圓心上的數字會往 LSB 為 0 的 MRN 捨入

問題回顧

了解上述內容後，讓我們回顧問題並釐清答案吧!

#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;
}

上述程式的輸出居然是 50002896，與實際答案 50005000 相差了 2104，雖然知道浮點數運算會有誤差，但 2104 是怎麼來的，你說的出來嗎?

透過 AI 撰寫以下程式分析該過程

/*
 * cumulative_error_print.c
 *
 * This program computes the cumulative sum from 1 to 10000 using both
 * float (simulating float32 arithmetic) and double (higher precision).
 * When the float and double cumulative results differ, the program prints
 * detailed information including:
 *   - Current number (num)
 *   - The current number as float (num_float)
 *   - Exact cumulative sum (exact_sum) computed using int64_t arithmetic
 *   - Cumulative sum computed with float (float32_sum)
 *   - LSB (least significant bit) of the float32 cumulative result
 *   - Current error (error = double_sum - float32_sum)
 *   - Change in error (error_diff compared to previous step)
 *
 * Compilation: gcc cumulative_error_print.c -lm
 * Execution: ./a.out > output.txt
 */

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

/* Use a union to access the bit-level representation of a float */
typedef union {
  float f;
  uint32_t u;
} FloatUnion;

/**
 * @brief Extract the least significant bit (LSB) of the mantissa from a float.
 *
 * @param value The float value.
 * @return The LSB of the 23-bit mantissa.
 */
unsigned int extract_lsb(float value) {
  FloatUnion fu;
  fu.f = value;
  /* Extract the mantissa (lower 23 bits) and then return its LSB */
  return (fu.u & 0x007FFFFF) & 1;
}

int main(void) {
  const int max_i = 10000; /* Loop from 1 to 10000 */
  int i;

  /* Exact cumulative sum computed with int64_t */
  int64_t exact_sum = 0;
  /* Cumulative sum using float (simulate float32 arithmetic) */
  float sum_f = 0.0f;
  /* Cumulative sum using double (higher precision, serves as theoretical value)
   */
  double sum_d = 0.0;

  double prev_error = 0.0;
  double current_error = 0.0;
  double error_diff = 0.0;

  /* Print header using tab-separated columns with left alignment */
  printf("%-8s\t%-10s\t%-12s\t%-12s\t%-4s\t%-8s\t%-8s\n", "num", "num_float",
         "exact_sum", "float32_sum", "LSB", "error", "error_diff");

  for (i = 1; i <= max_i; i++) {
    exact_sum += i;
    sum_d += i;               /* Accumulate using double (theoretical value) */
    sum_f = sum_f + (float)i; /* Accumulate using float (float32 simulation) */

    current_error = sum_d - (double)sum_f;

    /* Print information when the float and double results differ */
    if (fabs(current_error) > 0.0) {
      if (i == 1)
        error_diff = 0.0;
      else
        error_diff = current_error - prev_error;
      prev_error = current_error;

      /* Get the LSB of the float cumulative result */
      unsigned int lsb = extract_lsb(sum_f);
      printf("%-8d\t%-10.0f\t%-12" PRId64 "\t%-12.0f\t%-4d\t%-8.0f\t%-8.0f\n",
             i, (float)i, exact_sum, sum_f, lsb, current_error, error_diff);
    }
  }

  return 0;
}

結合上述得知僅當浮點數超過 $2^{24} 時，其 ULP 才會 > 1.0。
因此將 float32_sum 依照

2^{24}

(),

2^{25}

(33554432) 進行分界，得知 num 分別為 :

5793: 16776527 → 16782320
8192: 33549136 → 33557328

觀察這附近的輸出:

num     	num_float 	exact_sum   	float32_sum 	LSB 	error   	error_diff
5793    	5793      	16782321    	16782320    	0   	1       	1       
5794    	5794      	16788115    	16788114    	1   	1       	0       
5795    	5795      	16793910    	16793908    	0   	2       	1       
5796    	5796      	16799706    	16799704    	0   	2       	0       
5797    	5797      	16805503    	16805500    	0   	3       	1       
5798    	5798      	16811301    	16811298    	1   	3       	0       
5799    	5799      	16817100    	16817096    	0   	4       	1       
5800    	5800      	16822900    	16822896    	0   	4       	0
...
/* 以下是 2^25 範圍 */
8192    	8192      	33558528    	33557328    	0   	1200    	0       
8193    	8193      	33566721    	33565520    	0   	1201    	1       
8194    	8194      	33574915    	33573712    	0   	1203    	2       
8195    	8195      	33583110    	33581908    	1   	1202    	-1      
8196    	8196      	33591306    	33590104    	0   	1202    	0       
8197    	8197      	33599503    	33598300    	1   	1203    	1       
8198    	8198      	33607701    	33606496    	0   	1205    	2       
8199    	8199      	33615900    	33614696    	0   	1204    	-1      
8200    	8200      	33624100    	33622896    	0   	1204    	0       
8201    	8201      	33632301    	33631096    	0   	1205    	1       
8202    	8202      	33640503    	33639296    	0   	1207    	2       
8203    	8203      	33648706    	33647500    	1   	1206    	-1
.
.
.
9996    	9996      	49965006    	49962904    	0   	2102    	0       
9997    	9997      	49975003    	49972900    	1   	2103    	1       
9998    	9998      	49985001    	49982896    	0   	2105    	2       
9999    	9999      	49995000    	49992896    	0   	2104    	-1      
10000   	10000     	50005000    	50002896    	0   	2104    	0

在前半部分 (

[5793, 8192)

)，可以參考這張圖:

LSB 0 → LSB 1 的區域定義為 Region A (藍色處)
LSB 1 → LSB 0 的區域定義為 Region A (綠色處)

在藍色區域的中點(含)以左都會被捨入到左側的浮點數 (e.g., 步驟 2)，在綠色區域中點(含)以右都會被捨入到右側的浮點數。

誤差可以這樣統計:

對於任意等差數列和，可以找出其誤差差值 (error_diff) 循環規律，例如本題(遞增 1)，在
$[5793, 8192)$ 為 1, 0，循環長度為 ULP。在
$[8192, 11585)$ 為 0, 1, 2, -1，在
$[11585, 16384)$ 為 1, 2, 3, 4, -3, -2, -1, 0。每次循環的和是 ULP / 2。
統計 error_diff:
- $[5793, 8192)$ : 共 1199 組循環餘 1，單個循環 error_diff 和是 1，剩下的數字 8191 的 error 是 1，因此本階段共計 1200。
- $[8192, 10000)$ : 剛好 452 組，一組循環 2，共 904。
- 兩者和為 2104
區間內的等差數值數量不一定會整除 ULP，要計算最後一組剩下的 error_diff 加上前面每一組的 error_diff_sum

目前頂多是暴力解，有數值解嗎?算式如何列?

背後的原理

todo

由於浮點數間隔在一段距離內 (

2^{k} \to 2^{k + 1} - 1

) 是均勻分布的，又我們的目標數是遞增，所以

error 循環長度與 ULP 一致 (e.g., 0, 1, 2, -1，對應到 ULP = 4)

至此，我們知道累積的誤差來源正是因為捨入方向。IEEE 754 預設的捨入模式稱為銀行家捨入法 (Banker's rounding, 四捨六入五取偶)，詳見延伸閱讀: 捨入方法。

修正誤差 - Kahan summation algorithm

我們可以藉由 Kahan summation algorithm

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;
}

其原理是…

Kahan summation algorithm 的限制?

IEEE 754 Floating Point

浮點數標準

問題回顧

修正誤差 - Kahan summation algorithm

延伸閱讀: 捨入方法

延伸閱讀: Float16 V.S. bFloat16

延伸閱讀: 浮點數在 Linux Kernel 的使用

避免 overflow

Reference

IEEE 754 Floating Point

浮點數標準

問題回顧

修正誤差 - Kahan summation algorithm

延伸閱讀: 捨入方法

延伸閱讀: Float16 V.S. bFloat16

延伸閱讀: 浮點數在 Linux Kernel 的使用

避免 overflow

Reference

Read more

2025q1 Homework3 (kxo)

2025q1 Homework2 (quiz1+2)

2025q1 Homework1 (lab0)

ORB-SLAM3 for ROS2 on Ubuntu 24.04