程序计算精确圆周率π的方法

圆周率π = 3.14159265358… 是一个超越数，更是一个"正规数"——你的名字、生日、电话号码，甚至《西游记》整本书的编码，理论上都藏在π的无限小数位中。从古代的几何方法到现代计算机算法，人类计算π精度的能力在不断进步。直到2024年3月14日,美国加州计算机存储公司Solidigm宣布，他们耗时75天将圆周率已经计算到小数点后约105万亿位。

本文将带你从最简Python一行代码开始，逐层深入C语言高精度π计算的实现原理、代码优化技巧，并对比多种主流算法。文末给出可直接运行的"万进制加速版"代码，计算10000位仅需原版1/4的时间。

---

目录

1. 极简起步：Python + mpmath 一行搞定
2. Wallis级数推导与朴素C实现
3. 突破double瓶颈：GMP高精度库版本
4. 核心优化：数组模拟大数运算（10000位实战）
5. 性能翻倍：从10进制升级到万进制
6. 算法对比：Machin / Gauss-Legendre / BBP / Chudnovsky
7. 趣味扩展：自然底数e的高精度计算
8. 代码赏析：IOCCC获奖作品解读

---

1. 极简起步：Python + mpmath 一行搞定 {#1-极简起步}

如果你只想"用一下"π，Python 的 mpmath 库是最快的选择：

from mpmath import mp

mp.dps = 1000         # 设置精度为1000位
print(mp.pi)          # 输出1000位圆周率

mpmath 底层用 C 实现，内部使用的是 Chudnovsky 算法（后文会介绍），收敛极快。但如果想理解"高精度计算"的底层原理，我们从更基础的算法开始。

---

2. Wallis级数推导与朴素C实现 {#2-wallis级数}

2.1 公式推导

Wallis 乘积公式的级数展开形式为：

π/2 = 1 + 1/3 + (1/3)·(2/5) + (1/3)·(2/5)·(3/7) + (1/3)·(2/5)·(3/7)·(4/9) + ...

令每一项 zₙ 满足递推关系：

z₀ = 2
zₙ = zₙ₋₁ × a / b    （其中 a=1,2,3...， b=3,5,7...）

则：π = z₀ + z₁ + z₂ + …

2.2 朴素C实现（double精度，约14位）

#include <stdio.h>

int main()
{
    double x = 2.0, z = 2.0;
    int a = 1, b = 3;

    while (z > 1e-16)   // 【优化】原文用 z>0，理论上永远不会为0
    {
        z = z * a / b;
        x += z;
        a++;
        b += 2;
    }
    printf("π = %.14f\n", x);   // 输出: 3.14159265358979
    return 0;
}

局限：double 只有 8 字节（约 15-16 位有效数字），无法突破。

---

3. 突破 double 瓶颈：GMP 高精度库版本 {#3-gmp版本}

想用原生 C 获得更高精度，可以借助 GNU 的 GMP（大数库） + MPFR（高精度浮点库）：

#include <stdio.h>
#include <gmp.h>

int main()
{
    mpf_set_default_prec(512);  // 512比特 ≈ 154位十进制精度

    mpf_t x, z, term_a, term_b, eps;
    mpf_inits(x, z, term_a, term_b, eps, NULL);

    mpf_set_ui(x, 2);
    mpf_set_ui(z, 2);
    mpf_set_ui(term_a, 1);
    mpf_set_ui(term_b, 3);

    // 【优化】用 eps 控制精度，而非硬编码字符串比较
    mpf_set_str(eps, "1e-100", 10);

    int cnt = 0;
    while (mpf_cmp(z, eps) > 0)
    {
        mpf_mul(z, z, term_a);
        mpf_div(z, z, term_b);
        mpf_add(x, x, z);
        mpf_add_ui(term_a, term_a, 1);
        mpf_add_ui(term_b, term_b, 2);
        cnt++;
    }

    gmp_printf("迭代 %d 次\nπ = %.100Ff\n", cnt, x);

    mpf_clears(x, z, term_a, term_b, eps, NULL);
    return 0;
}
// 编译: gcc main.c -lgmp -lmpfr

但 GMP 终究是"外挂"，下面我们来手写一个纯C的高精度方案。

---

4. 核心优化：数组模拟大数运算（10000位实战） {#4-数组模拟大数}

4.1 基本原理

用一个 char 数组，每个元素存储 0-9 的一位十进制数，模拟手算竖式的 乘法进位 和 除法试商 过程：

数组布局: [x₀][x₁].[x₂][x₃]...[xₙ]
           ↑ 整数部分  ↑ 小数点    ↑ 小数部分

4.2 优化版代码（含详细注释）

#include <stdio.h>
#include <string.h>

int main()
{
    const int ARRSIZE = 10100;  // 数组大小（多留一些防止溢出）
    const int DISPCNT = 10000;  // 实际显示位数

    // x[0]备用, x[1]为整数部分, x[2..n]为小数部分
    char x[10100], z[10100];
    int a = 1, b = 3;
    int i, run_flag, cnt = 0;

    memset(x, 0, sizeof(x));
    memset(z, 0, sizeof(z));
    x[1] = 2;
    z[1] = 2;

    // 【优化】使用 do-while 避免每次多算一轮
    do
    {
        // --- 步骤1: z = z * a ---
        int carry = 0;
        for (i = ARRSIZE - 1; i >= 0; i--)
        {
            int tmp = z[i] * a + carry;
            z[i] = tmp % 10;
            carry = tmp / 10;
        }

        // --- 步骤2: z = z / b ---
        carry = 0;
        for (i = 0; i < ARRSIZE; i++)
        {
            int tmp = z[i] + carry * 10;
            z[i] = tmp / b;
            carry = tmp % b;
        }

        // --- 步骤3: x = x + z ---
        // 同时检测 z 是否全为0（收敛判断）
        run_flag = 0;
        for (i = ARRSIZE - 1; i > 0; i--)
        {
            int sum = x[i] + z[i];
            x[i] = sum % 10;
            x[i - 1] += sum / 10;
            run_flag |= z[i];   // 只要有非零位就继续
        }

        a++;
        b += 2;
        cnt++;
    } while (run_flag && cnt < 100000);

    // --- 输出 ---
    printf("迭代 %d 次\n", cnt);
    printf("π = %d.", x[1]);
    for (i = 0; i < DISPCNT; i++)
    {
        if (i > 0 && i % 100 == 0)  printf("\n");
        printf("%d", x[i + 2]);
    }
    printf("\n");
    return 0;
}

运行结果：迭代约 33538 次 即可达到 10000 位精度。该算法约 3.3 轮迭代产生一位正确数字。

---

5. 性能翻倍：从10进制升级到万进制 {#5-万进制加速}

5.1 为什么能提速？

上述代码每个数组元素只存 0-9，每次迭代要对 ~10000 个元素各做一次乘除——相当于 10000 次单精度运算。

如果把每个元素改为存 0-9999（万进制），数组长度缩为原来的 1/4，速度直接提升约 4倍：

#include <stdio.h>
#include <string.h>

#define RADIX  10000    // 万进制基数
#define WIDTH  4        // 每个元素对应4位十进制

int main()
{
    const int ARRSIZE = 2600;   // 万进制数组（2600×4 ≈ 10400位）
    const int DISPCNT  = 10000;

    int x[ARRSIZE], z[ARRSIZE];
    int a = 1, b = 3;
    int cnt = 0, run_flag;

    memset(x, 0, sizeof(x));
    memset(z, 0, sizeof(z));
    x[1] = 2;
    z[1] = 2;

    do
    {
        // --- z = z * a / b ---
        long long carry = 0;
        int i;
        for (i = ARRSIZE - 1; i >= 0; i--)
        {
            carry += (long long)z[i] * a;
            z[i] = (int)(carry % RADIX);
            carry /= RADIX;
        }

        carry = 0;
        for (i = 0; i < ARRSIZE; i++)
        {
            carry = carry * RADIX + z[i];
            z[i] = (int)(carry / b);
            carry %= b;
        }

        // --- x = x + z ---
        run_flag = 0;
        carry = 0;
        for (i = ARRSIZE - 1; i >= 0; i--)
        {
            carry += (long long)x[i] + z[i];
            x[i] = (int)(carry % RADIX);
            carry /= RADIX;
            run_flag |= z[i];
        }

        a++;
        b += 2;
        cnt++;
    } while (run_flag && cnt < 100000);

    // --- 格式化输出 ---
    printf("迭代 %d 次（万进制加速版）\n", cnt);
    printf("π = %d.", x[1]);
    for (int i = 0; i < DISPCNT; i++)
    {
        if (i > 0 && i % 100 == 0) printf("\n");
        // 从万进制数组中取出每一位十进制数字
        int idx = i / WIDTH + 2;
        int pos = WIDTH - 1 - (i % WIDTH);
        int digit = x[idx];
        while (pos--) digit /= 10;
        printf("%d", digit % 10);
    }
    printf("\n");
    return 0;
}

性能对比（计算10000位）：

版本	数组大小	每次迭代运算量	相对耗时
原版（十进制）	10100	~10100次单精度运算	1.0×
万进制加速版	2600	~2600次32位运算	~0.25×

---

6. 算法对比：哪种算法最快？ {#6-算法对比}

Wallis 级数虽然实现简单，但收敛速度并不理想（约 3.3 轮/位）。以下是主流π计算算法的对比：

算法	每次迭代增加位数	10000位迭代次数	复杂度	实现难度
Wallis 级数（本文）	~0.3 位	~33,000	O(n²)	★☆☆☆☆
Machin 公式 (1706)	~1.4 位	~7,100	O(n²)	★★☆☆☆
Gauss-Legendre (1975)	翻倍	~14	O(n log²n)	★★★☆☆
BBP 公式 (1995)	~1.2 位	~8,300	O(n²)	★★★☆☆
Chudnovsky (1988)	~14 位	~715	O(n log³n)	★★★★★

6.1 Machin 公式（推荐入门优化）

π/4 = 4·arctan(1/5) - arctan(1/239)

收敛速度约 1.4 位/项，比 Wallis 快 4 倍以上，且实现难度适中：

// arctan(1/x) 的级数展开：1/x - 1/(3x³) + 1/(5x⁵) - ...
// 套用数组大数运算即可，留作读者练习

6.2 Gauss-Legendre 算法（二次收敛，王者速度）

该算法每次迭代精度翻倍，计算 10000 位仅需 14 次迭代：

初始化: a₀ = 1,  b₀ = 1/√2,  t₀ = 1/4,  p₀ = 1
迭代:
  aₙ₊₁ = (aₙ + bₙ) / 2
  bₙ₊₁ = √(aₙ × bₙ)
  tₙ₊₁ = tₙ - pₙ × (aₙ - aₙ₊₁)²
  pₙ₊₁ = 2 × pₙ
结果: π ≈ (aₙ + bₙ)² / (4 × tₙ)

这是目前不使用FFT乘法前提下最快的π算法，强烈推荐作为进阶学习目标。

6.3 Chudnovsky 算法（工业级）

mpmath、Mathematica 等专业软件的首选，每项增加约 14 位精度：

1/π = 12 × Σ (-1)ᵏ × (6k)! × (13591409 + 545140134k)
           ─────────────────────────────────────────
           (3k)! × (k!)³ × 640320³ᵏ⁺³/²

需要配合**二分拆分（Binary Splitting）**技术才能高效实现。

---

7. 趣味扩展：自然底数 e 的高精度计算 {#7-扩展}

同样的数组模拟思路也适用于 e = 2.71828…：

#include <stdio.h>

int main()
{
    const int TIMES = 1000;
    int f[TIMES + 1];
    int n, d, e_digit = 0;

    for (n = TIMES; n > 0; n--)
        f[n] = 10;

    printf("e = 2.");
    for (int i = 0; i < TIMES; i++)
    {
        d = 0;
        for (n = TIMES; n > 0; n--)
        {
            d = d + f[n] * 10;
            f[n] = d % n;
            d = d / n;
        }
        printf("%d", (e_digit + d / 10) % 10);
        e_digit = d % 10;
    }
    printf("\n");
    return 0;
}

原理基于：e = 1 + 1/1! + 1/2! + 1/3! + ...，利用秦九韶算法（Horner’s method）变形实现。

---

8. 代码赏析：IOCCC 获奖作品解读 {#8-代码赏析}

8.1 800位圆周率（经典一行流）

某年 IOCCC（国际C语言代码混淆竞赛） 的获奖作品，仅 3行代码 输出800位π：

#include <stdio.h>
long a = 10000, b, c = 2800, d, e, f[2801], g;
int 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)
            ;
}

这段代码实现的是 Spigot 算法（逐位输出），利用"整数只留余数"的方式避免了高精度浮点运算。

8.2 Bellard 的梅森素数计算器

Bellard 在 IOCCC 2000 上的获奖作品，使用 NTT（数论变换） 在极短时间内验证了第38个梅森素数（2⁶⁹⁷²⁵⁹³−1，有 2,098,960 位）。这展示了"用整数运算模拟大数"能达到的惊人高度。

---

一句话：理解原理用 Wallis 级数 + 数组大数；追求速度用 Gauss-Legendre；工程落地直接调 GMP/mpmath。

---

参考：

• Wallis product - Wikipedia
• Chudnovsky algorithm - Wikipedia
• Gauss-Legendre algorithm - Wikipedia
• IOCCC 官方网站