C語言高精度計算：圓周率前 1000 位數值生成器

馬赫恩類公式 + 大數運算：純 C 實作圓周率任意位數計算

使用高斯/史托默公式結合大整數陣列模擬，計算圓周率任意位數。這是一個 1999 年的經典作品，展示純整數運算實現高精度數學常數計算。

核心數學公式

高斯公式 (預設)：

π = 48·arctan(1/18) + 32·arctan(1/57) - 20·arctan(1/239)

史托默公式 (驗算用)：

π = 24·arctan(1/8) + 8·arctan(1/57) + 4·arctan(1/239)

利用反正切級數展開：arctan(1/x) = 1/x - 1/(3x³) + 1/(5x⁵) - ...

程式架構解析

// 核心資料結構：長整數陣列模擬大數
// 每個元素儲存 1 位十進制數字 (0-9)，低索引 = 低位
long *pi, *t;   // pi: 累積結果, t: 暫存單項計算
int d;          // 目標位數 + 5 守衛位

// 主流程：
// 1. 解析命令列參數 (位數、公式選擇)
// 2. 分配陣列 (calloc 歸零)
// 3. 對每個 arctan 項：
//    a. 初始化 t = 係數 / 分母
//    b. 長除法迴疊計算級數各項
//    c. 正負交替累加到 pi
// 4. 處理進位/借位 (修正每位為 0-9)
// 5. 格式化輸出 (每 5 位分組，每 65 位換行)

關鍵演算法細節

長除法模擬 (核心)

// 計算 t = t / n (大數除以小數)
int r = 0;
for (int i = highest_digit; i >= 0; i--) {
    r = r * 10 + t[i];      // 帶下下一位
    t[i] = r / n;           // 商存回
    r = r % n;              // 餘數留給下一位
}

級數項遞推關係

arctan(1/x) = Σ (-1)^k / ((2k+1) * x^(2k+1))
第 k 項 = 前一項 / (-x²) / (2k+1) * (2k-1)
實作中：分母平方後，每兩步除以 (2k+1) 並交替加減

執行範例

$ ./pi 1000
Nature (R) Pi value compute Program (C) Tue 1999.11.30
pi= 48 * arctg(1/18) + 32 * arctg(1/57) - 20 * arctg(1/239) [Gauss]
pi= 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 ...
DIGITS: 1000

$ ./pi 5000 Stomer
pi= 24 * arctg(1/8) + 8 * arctg(1/57) + 4 * arctg(1/239) [Stomer]
... (5000 位輸出)
DIGITS: 5000

效能與精度特性

現代化改良方向

• 基底優化：改用 10⁹ 為基底，每元素存 9 位，陣列縮小 9 倍
• FFT 乘法：大數乘法用 FFT 從 O(n²) 降為 O(n log n)
• Chudnovsky 公式：收斂極快，每項增加 ~14 位
• 平行化：多執行緒計算不同 arctan 項
• 記憶體對齊：使用 aligned_alloc 配合 SIMD

歷史意義

這類「手工大數運算」程式展示了在無任意精度庫 (GMP, MPFR) 時代，如何用基礎 C 語言實現高精度計算。現代生產環境建議直接使用 GMP (mpz_t, mpf_t) 或 MPFR，效能、正確性、維護性皆優於手寫實作。