1989 - Kahan 的 IEEE 754（浮點數的統一標準）#

import numpy as np
import struct

# 1) 浮點數不能精確表示（二進位的悲劇）
print("浮點數的表示（二進位）：")
print(f"  0.1 + 0.2 = {0.1 + 0.2!r}")
print(f"  0.1 + 0.2 == 0.3？{0.1 + 0.2 == 0.3}")
print(f"  0.1 + 0.2 - 0.3 = {0.1 + 0.2 - 0.3!r}（誤差殘留）\n")

# 2) IEEE 754 的位元結構（32 位元單精度）
def float_bits(x):
    """顯示 IEEE 754 的位元（32 位元）"""
    [b] = struct.unpack('>I', struct.pack('>f', x))
    sign = (b >> 31) & 1
    exp = (b >> 23) & 0xFF
    frac = b & 0x7FFFFF
    return sign, exp, frac

print("IEEE 754 位元結構（32 位元單精度）：")
for x in [1.0, 2.0, 0.5, -1.0, 0.1]:
    s, e, f = float_bits(x)
    print(f"  {x}: 符號 = {s}, 指數 = {e} (127+{e-127}), 尾數 = {f:023b}")
print()

# 3) 特殊值：∞、NaN、非正規數
print("特殊值（Kahan 的設計）：")
inf = np.float64('inf'); nan = np.float64('nan')
print(f"  1.0 / 0.0 = {np.float64(1.0) / np.float64(0.0)}（+∞）")
print(f"  -1.0 / 0.0 = {np.float64(-1.0) / np.float64(0.0)}（−∞）")
print(f"  0.0 / 0.0 = {np.float64(0.0) / np.float64(0.0)}（NaN）")
print(f"  ∞ - ∞ = {inf - inf}（NaN）")
print(f"  NaN == NaN？{nan == nan}（NaN 不等於自己！）")
# 非正規數
smallest_normal = np.finfo(np.float64).tiny
smallest_subnormal = np.nextafter(0.0, 1.0)
print(f"  最小正規數 = {smallest_normal!r}")
print(f"  最小非正規數 = {smallest_subnormal!r}（漸進下溢）\n")

# 4) 最近偶數捨入（Kahan 的關鍵設計）
print("最近偶數捨入（round-to-nearest-even）：")
# 2.5 應捨入到 2（偶數），不是 3
import decimal
from decimal import Decimal, ROUND_HALF_EVEN
for x in ['1.5', '2.5', '3.5', '0.5']:
    d = Decimal(x).quantize(Decimal('1'), rounding=ROUND_HALF_EVEN)
    print(f"  {x} → {d}（捨入到偶數）")
print("  → 2.5 → 2（不是 3！）：統計無偏（誤差不累積）\n")

# 5) 捨入誤差的累積（災難的模擬）
print("捨入誤差的累積（Kahan 總和法 vs 天真總和）：")
def naive_sum(arr):
    """天真總和（誤差累積）"""
    s = 0.0
    for x in arr:
        s += x
    return s

def kahan_sum(arr):
    """Kahan 總和法（補償求和，誤差極小）"""
    s = 0.0
    c = 0.0   # 補償項
    for x in arr:
        y = x - c
        t = s + y
        c = (t - s) - y   # 捕捉丟失的位元
        s = t
    return s

# 大量小數的總和
n = 10**7
arr = [0.1] * n   # 0.1 不能精確表示
naive = naive_sum(arr)
kahan = kahan_sum(arr)
print(f"  {n} 個 0.1 的總和：")
print(f"  天真總和 = {naive!r}")
print(f"  Kahan 總和 = {kahan!r}")
print(f"  理論值 = {n*0.1!r}")
print(f"  天真誤差 = {abs(naive - n*0.1):.2e}")
print(f"  Kahan 誤差 = {abs(kahan - n*0.1):.2e}")
print("  → Kahan 總和法把誤差降到最低（1989 圖靈獎的方法）")

浮點數的表示（二進位）：
  0.1 + 0.2 = 0.30000000000000004
  0.1 + 0.2 == 0.3？False
  0.1 + 0.2 - 0.3 = 5.551115123125783e-17（誤差殘留）

IEEE 754 位元結構（32 位元單精度）：
  1.0: 符號 = 0, 指數 = 127 (127+0), 尾數 = 00000000000000000000000
  2.0: 符號 = 0, 指數 = 128 (127+1), 尾數 = 00000000000000000000000
  0.5: 符號 = 0, 指數 = 126 (127+-1), 尾數 = 00000000000000000000000
  -1.0: 符號 = 1, 指數 = 127 (127+0), 尾數 = 00000000000000000000000
  0.1: 符號 = 0, 指數 = 123 (127+-4), 尾數 = 10011001100110011001101

特殊值（Kahan 的設計）：
  1.0 / 0.0 = inf（+∞）
  -1.0 / 0.0 = -inf（−∞）
  0.0 / 0.0 = nan（NaN）
  ∞ - ∞ = nan（NaN）
  NaN == NaN？False（NaN 不等於自己！）
  最小正規數 = 2.2250738585072014e-308
  最小非正規數 = 5e-324（漸進下溢）

最近偶數捨入（round-to-nearest-even）：
  1.5 → 2（捨入到偶數）
  2.5 → 2（捨入到偶數）
  3.5 → 4（捨入到偶數）
  0.5 → 0（捨入到偶數）
  → 2.5 → 2（不是 3！）：統計無偏（誤差不累積）

捨入誤差的累積（Kahan 總和法 vs 天真總和）：
  10000000 個 0.1 的總和：
  天真總和 = 999999.9998389754
  Kahan 總和 = 1000000.0
  理論值 = 1000000.0
  天真誤差 = 1.61e-04
  Kahan 誤差 = 0.00e+00
  → Kahan 總和法把誤差降到最低（1989 圖靈獎的方法）