C語言浮點數精度問題深度解析：Fortran 遷移、自定義浮點型別與高精度計算策略

問題背景

將 Fortran 程式移植到 C，原程式在第 9 位小數截斷，C 語言用 pow 計算冪次出現第 14 位小數，第 8 位即產生進位誤差。程式被呼叫上萬次，要求精度「絲毫不差」。

核心原因：IEEE 754 二進位浮點數本質

• 十進位小數多數無法精確用二進位表示：0.1、0.01、0.001 在二進位皆為無限循環小數。
• float (單精度) 約 6-7 位十進位有效數字，double (雙精度) 約 15-16 位。
• pow、exp、log 等超越函數由標準庫實作，誤差通常在 0.5-1 ULP（最後一位單位），累積上萬次誤差不可控。
• Fortran 某些編譯器支援十進位浮點或定點運算，C 標準庫不保證。

解決方案對比

推薦方案一：定點整數運算（最工程化、最高效）

若業務需求為「固定小數位（如 9 位）精確計算」，直接用 int64_t 儲存「億分之一」單位：

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

#define SCALE 1000000000L  // 10^9，對應 9 位小數

// 定點加法
int64_t fixed_add(int64_t a, int64_t b) { return a + b; }

// 定點減法
int64_t fixed_sub(int64_t a, int64_t b) { return a - b; }

// 定點乘法：結果需除以 SCALE，注意溢位與四捨五入
int64_t fixed_mul(int64_t a, int64_t b) {
    __int128 tmp = (__int128)a * b;  // GCC/Clang 支援 128 位中間值
    return (int64_t)((tmp + SCALE/2) / SCALE);  // 四捨五入
}

// 定點除法：被除數先乘 SCALE
int64_t fixed_div(int64_t a, int64_t b) {
    __int128 tmp = (__int128)a * SCALE;
    return (int64_t)((tmp + b/2) / b);  // 四捨五入
}

// 從 double 轉換（僅初始化用）
int64_t from_double(double x) { return (int64_t)(x * SCALE + 0.5); }

// 輸出顯示
void print_fixed(int64_t x) {
    int64_t int_part = x / SCALE;
    int64_t frac_part = x % SCALE;
    if (frac_part < 0) frac_part = -frac_part;
    printf("%" PRId64 ".%09" PRId64 "\n", int_part, frac_part);
}

推薦方案二：十進位浮點 _Decimal64 (GCC/Clang 擴充，C23 標準化)

// 需 GCC/Clang，編譯加 -std=c23 或 -fdecimal-float
#include <stdio.h>

_Decimal64 a = 0.1dd;   // 十進位字面常數後綴 dd
_Decimal64 b = 0.2dd;
_Decimal64 c = a + b;   // 精確為 0.3，無二進位誤差
printf("%.20Df\n", c);  // 輸出 0.30000000000000000000

推薦方案三：GMP/MPFR 任意精度（最高精度、最重）

// 需連結 -lgmp -lmpfr
#include <mpfr.h>

mpfr_t x, y, z;
mpfr_init2(x, 256);  // 256 位精度
mpfr_init2(y, 256);
mpfr_init2(z, 256);
mpfr_set_str(x, "0.1", 10, MPFR_RNDN);
mpfr_set_str(y, "0.2", 10, MPFR_RNDN);
mpfr_add(z, x, y, MPFR_RNDN);  // 精確十進位加法
mpfr_printf("%.50Rf\n", z);    // 0.30000000000000000000000000000000000000000000000000

原回答「自定義結構」評析

typedef struct MyFloat {
    long a;  // 尾數
    long n;  // 10 的冪次
} MyFloat;

• 優點：概念清晰，尾數為整數運算精確，冪次分離避免大數溢位。
• 缺點：需自行實作加減乘除、比較、輸入輸出、正規化（尾數範圍調整）、捨入模式；long 範圍有限（通常 64 位），大數仍會溢位；效能遠遜原生浮點。
• 結論：適合教學理解科學記數法、或極簡場景；工程專案建議用成熟庫。

Fortran 遷移具體建議

1. 確認精度需求：業務真正需要「位元級完全一致」還是「業務邏輯正確」？前者極難，後者可接受誤差界限。
2. 統一精度策略：全系統統一用 double (雙精度) 或 long double (擴展精度)，避免混用 float。
3. 替換 pow：整數冪次用遞推乘法；小數冪次若精度敏感，考慮查表或多項式逼近（如 exp(y*log(x))）。
4. 累積誤差補償：大量加法用 Kahan 求和算法：

double kahan_sum(const double *arr, size_t n) {
    double sum = 0.0, c = 0.0;
    for (size_t i = 0; i < n; ++i) {
        double y = arr[i] - c;
        double t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }
    return sum;
}

總結

C 語言原生浮點無法保證十進位精度「絲毫不差」。工程上三大主流方案：

1. 定點整數（固定小數位、高效能、絕對精確）——首選。
2. 十進位浮點 _Decimal64/128（標準化、語意正確、硬體支援中）。
3. GMP/MPFR（任意精度、科學計算基準、重型依賴）。

自定義結構可作教學或極輕量場景，生產環境請用成熟庫。