C 語言計算圓周率 π 的極致壓縮程式：Spigot 算法完整解析

編輯：關於C語言

前言：只用幾行代碼算出上萬位 π

這段傳奇的 C 代碼僅 5 行、不足 200 字元，卻能一次性輸出圓周率 π 的 2400+ 位數字。它出自 1990 年代國際混亂 C 程式碼大賽 (IOCCC)，展示了 Spigot 算法（噴嘴算法） 以整數運算逐位產生 π 數字的驚人威力。本文將完整還原代碼、逐行解析算法原理、記憶體佈局，並提供可讀性重構版本。

原始極致壓縮代碼

long a=10000,b,c=2800,d,e,f[2801],g;
main() {
    for(;b-c;) f[b++]=a/5;
    for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
        for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

格式化後的可讀版本

#include 

#define BASE 10000      // 每次輸出 4 位十進制數字
#define DIGITS 2800     // 目標位數（約 2800 * log10(10000) ≈ 11200 位）
#define ARRAY_SIZE 2801 // 工作陣列大小

long a = BASE;          // 基數 10000
long b, c = DIGITS;     // 迴圈計數器
long d, e;              // 暫存變數
long f[ARRAY_SIZE];     // 工作陣列：儲存分數係數
long g;                 // 分母迴圈變數

int main() {
    // === 初始化階段 ===
    // f[i] = a/5 = 2000，對應級數展開的初始係數
    for (; b < c; ) {
        f[b++] = a / 5;  // 2000
    }

    // === 主計算迴圈 ===
    // 每次迴圈產生 4 位數字 (c -= 14 對應 14 項級數項)
    for (; d = 0, g = c * 2; c -= 14, printf("%.4d", e + d / a), e = d % a) {
        // 內層迴圈：從高項到低項進行「進位傳播」
        for (b = c; d += f[b] * a, f[b] = d % --g, d /= g--, --b; d *= b);
    }
    return 0;
}

算法核心原理：Spigot 算法

數學基礎：π 的級數展開

算法基於以下級數（源自阿基米德/萊布尼茨類公式的變體）：

$$ \pi = \sum_{k=0}^{\infty} \frac{2^k (k!)^2}{(2k+1)!} = 2 + \frac{1}{3}\left(2 + \frac{2}{5}\left(2 + \frac{3}{7}\left(2 + \cdots\right)\right)\right) $$

這是一個連分數形式，適合從高項向低項遞推計算。

整數運算模擬分數運算

核心洞見：用整數陣列 f[] 存儲各項的「分子係數」，基數 a=10000 作為縮放因子，將所有分數運算轉為整數乘除：

f[b] * a：將係數放大基數倍
d += ...：累加進位
f[b] = d % --g：計算當前項新係數（餘數）
d /= g--：進位傳給下一項（商）
d *= b：乘以項索引（對應階乘關係）

記憶體佈局與變數角色對照表

變數	角色	初值/範圍
`a`	輸出基數 (10^4)	10000
`c`	剩餘項數計數器	2800 → 0
`f[]`	工作陣列，存各項係數	初始全 2000
`g`	分母迴圈變數 (2c, 2c-1...)	動態變化
`d`	進位累加器	每輪歸零
`e`	上一輪輸出的餘數	跨輪傳遞

執行流程逐輪追蹤（前 3 輪）

初始: f[0..2799] = 2000, c=2800, e=0

第 1 輪 (c=2800, g=5600):
  內層從 b=2799 遞減至 0:
    d 累加進位, f[b] 更新為餘數, d 為商傳遞
  內層結束: d ≈ 31415926..., e = d % 10000
  輸出: printf("%.4d", e + d/a) → 3141
  c=2786, e 保留供下輪

第 2 輪 (c=2786, g=5572):
  類似計算...
  輸出: 5926

第 3 輪:
  輸出: 5358
...

輸出結果驗證（前 100 位）

3.14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679
...

與標準 π 值完全吻合。完整程式可輸出約 11200 位數字（受限於 long 32 位溢出邊界）。

現代化改良版本（支援任意精度）

#include 
#include 
#include   // GNU Multiple Precision Arithmetic Library

void pi_spigot_mpz(int digits) {
    int base = 10000;
    int terms = digits * 14 / 4 + 10;  // 經驗公式
    mpz_t *f = calloc(terms, sizeof(mpz_t));
    mpz_t a, d, e, g, tmp;
    mpz_init_set_ui(a, base);
    mpz_init(d); mpz_init(e); mpz_init(g); mpz_init(tmp);

    for (int i = 0; i < terms; i++) mpz_init_set_ui(f[i], base / 5);

    for (int c = terms; c > 0; c -= 14) {
        mpz_set_ui(d, 0);
        mpz_set_ui(g, 2 * c);
        for (int b = c - 1; b >= 0; b--) {
            mpz_mul_ui(tmp, f[b], base);
            mpz_add(d, d, tmp);
            mpz_sub_ui(g, g, 1);
            mpz_mod(f[b], d, g);
            mpz_div(d, d, g);
            mpz_mul_ui(d, d, b);
        }
        mpz_div(tmp, d, a);
        mpz_add(tmp, tmp, e);
        gmp_printf("%04Zd", tmp);
        mpz_mod(e, d, a);
    }
    // 清理...
}