C 語言計算圓周率：經典反正切級數實現與多公式支援

程式由來與特色

這是一個 1999 年編寫的 C 語言圓周率計算程式，純數學運算，無平台相依性，可在 Linux、Unix、Windows 下編譯執行。程式支援三種經典公式：

• Gauss（預設）：π = 48·arctan(1/18) + 32·arctan(1/57) - 20·arctan(1/239)
• Störmer（第二參數非零）：π = 24·arctan(1/8) + 8·arctan(1/57) + 4·arctan(1/239)
• 命令列參數控制位數與公式選擇

核心算法解析

泰勒級數與長除法

反正切級數：

arctan(1/d) = 1/d - 1/(3d³) + 1/(5d⁵) - 1/(7d⁷) + ...

程式以 long 陣列模擬高精度整數，每元素存一位十進制數（基數 10），低位在低索引。核心運算為長除法：

for (i = j; i >= 0; i--) {
    r = (m = 10 * r + t[i]) % n;
    t[i] = m / n;
}

其中 t 為被除數，n 為除數，r 為餘數。商直接累加/減到結果陣列 pi。

雙重迴圈結構

外層迴圈遍歷公式三項（p = 0..2）。內層第一階段計算首項 c/d；第二階段 while (j > 0) 迭代計算後續項：先平方分母 n *= n，再長除得 1/(d²ᵏ⁺¹)，最後除以奇數步長 s += 2，正負交替累加。

進位修正與輸出

所有項累加完成後，pi 陣列可能包含負數或 ≥10 的值。最終統一進位修正：

for (i = 0; i <= d; pi[i++] = r) {
    n = (m = pi[i] + n) / 10;
    if ((r = m % 10) < 0) r += 10, n--;
}

輸出格式：每行 65 位，每 5 位一組，最後顯示總位數。

執行範例

> pi.exe 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

歷史註記與現代視角

• 原程式假設 int 為 16 位，現代編譯器多為 32 位，但邏輯不受影響。
• 基數 10 效率較低，現代高精度庫多採用基數 2³² 或 10⁹ 配合 FFT 乘法。
• 程式展示了「以空間換時間」的陣列模擬大整數思想，是理解任意精度運算的優秀教學案例。