高斯-勒让德算法 （Gauss-Legendre Algorithm），也称为Brent-Salamin 算法 ，是一种利用算术-几何平均（Arithmetic-Geometric Mean, AGM）快速计算圆周率 π π 二次收敛 著称——每次迭代使有效数字翻倍，仅需 25 次迭代即可计算出超过 4500 万位 π π
圆周率 π π π π
高斯-勒让德算法代表了 π π 一次革命性飞跃 ——它利用算术-几何平均（AGM）的二次收敛特性，用极少的迭代次数即可获得极高精度的 π π
特性
说明
收敛速度 二次收敛（quadratic convergence），每步有效数字翻倍
迭代次数 25 次可达 4500 万位精度
核心数学 算术-几何平均（AGM）、椭圆积分
开发者 高斯（Carl Friedrich Gauss）、勒让德（Adrien-Marie Legendre）；1976 年由 Richard Brent 和 Eugene Salamin 独立重新发现
历史意义 1976 年首次突破百万位 π π
1799 年，22 岁的高斯在研究椭圆积分时，偶然发现了一个惊人的现象：对于任意两个正数 a a b b
高斯在他的日记中写道：
"算术-几何平均的研究揭示了椭圆积分与初等函数之间深刻的关系，这是我最重要的发现之一。"
与此同时，法国数学家勒让德也在独立研究椭圆积分，发现了一类特殊的积分关系。后来历史学家证实，勒让德实际上比高斯更早描述了 AGM 与 π π 高斯-勒让德算法 。
1976 年，澳大利亚数学家 Richard Brent 和美国数学家 Eugene Salamin 分别独立地发现，利用 AGM 可以构造一个二次收敛 的 π π
Brent, R. P. (1976). Multiple-precision zero-finding methods and the complexity of elementary function evaluation
Salamin, E. (1976). Computation of \pi using arithmetic-geometric mean
这一发现使 π π 数千位 一跃突破百万位 。
对于两个正实数 a a b b a ≥ b > 0 a ≥ b > 0
a n + 1 = a n + b n 2 （算术平均） a n + 1 = 2 a n + b n （算术平均）
b n + 1 = a n b n （几何平均） b n + 1 = a n b n （几何平均）
初始 a 0 = a a 0 = a b 0 = b b 0 = b AGM ( a , b ) AGM ( a , b )
让我们用一个简单例子来理解 AGM 的收敛过程。
取 a 0 = 10 a 0 = 10 b 0 = 1 b 0 = 1
迭代次数 n n
a n a n b n b n a n − b n a n − b n
0
10.000000000000000
1.000000000000000
9.000000000000000
1
5.500000000000000
3.162277660168380
2.337722339831620
2
4.331138830084190
4.169652381589939
0.161486448494251
3
4.250395605837064
4.249934869380331
0.000460736456733
4
4.250165237608698
4.250165237259122
3.5 × 10 − 10 3.5 × 1 0 − 10
5
4.250165237433910
4.250165237433910
< 10 − 15 < 1 0 − 15
观察可知：
第 0 步：两个数相差 9，差异巨大
第 3 步：差异缩小到 4.6 × 10 − 4 4.6 × 1 0 − 4
第 4 步：差异缩小到 3.5 × 10 − 10 3.5 × 1 0 − 10
第 5 步：两数在计算机双精度范围内完全相等
关键观察 ：AGM 的收敛速度远快于单独计算算术平均或几何平均，这是二次收敛 的体现——每步的有效数字精准地翻倍。
高斯证明了 AGM 与第一类完全椭圆积分 K ( k ) K ( k )
K ( k ) = ∫ 0 π / 2 d θ 1 − k 2 sin 2 θ = π 2 ⋅ AGM ( 1 , 1 − k 2 ) K ( k ) = ∫ 0 π /2 1 − k 2 sin 2 θ d θ = 2 ⋅ AGM ( 1 , 1 − k 2 ) π
这一关系正是高斯-勒让德算法的数学根基：通过 AGM 来高效计算椭圆积分，进而得到 π π 。
初始化：
a 0 = 1 , b 0 = 1 2 , t 0 = 1 4 , p 0 = 1 a 0 = 1 , b 0 = 2 1 , t 0 = 4 1 , p 0 = 1
迭代（对 n ≥ 0 n ≥ 0
a n + 1 = a n + b n 2 a n + 1 = 2 a n + b n
b n + 1 = a n b n b n + 1 = a n b n
t n + 1 = t n − p n ( a n − a n + 1 ) 2 t n + 1 = t n − p n ( a n − a n + 1 ) 2
p n + 1 = 2 p n p n + 1 = 2 p n
π π
π n ≈ ( a n + 1 + b n + 1 ) 2 4 t n + 1 π n ≈ 4 t n + 1 ( a n + 1 + b n + 1 ) 2
该算法的正确性源于勒让德发现的椭圆积分恒等式：
π = [ AGM ( a 0 , b 0 ) ] 2 a 0 2 − ∑ n = 0 ∞ 2 n ( a n − a n + 1 ) 2 π = a 0 2 − ∑ n = 0 ∞ 2 n ( a n − a n + 1 ) 2 [ AGM ( a 0 , b 0 ) ] 2
其中 a 0 = 1 a 0 = 1 b 0 = 1 / 2 b 0 = 1/ 2 t n t n p n p n 2 n 2 n
让我们用实际数值展示每一步的计算过程。
初始值：
a 0 = 1.0000000000000000 a 0 = 1.0000000000000000
b 0 = 1 2 ≈ 0.7071067811865475 b 0 = 2 1 ≈ 0.7071067811865475
t 0 = 0.2500000000000000 t 0 = 0.2500000000000000
p 0 = 1 p 0 = 1
初始近似：
π 0 ≈ ( 1 + 0.7071067811865475 ) 2 4 × 0.25 = 2.914213562373095 1 = 2.914213562373095 π 0 ≈ 4 × 0.25 ( 1 + 0.7071067811865475 ) 2 = 1 2.914213562373095 = 2.914213562373095
误差约 0.227（误差较大，因为仅用初始值的近似还很粗糙）。
第 1 次迭代：
a 1 = 1 + 0.7071067811865475 2 = 0.8535533905932737 a 1 = 2 1 + 0.7071067811865475 = 0.8535533905932737
b 1 = 1 × 0.7071067811865475 = 0.7071067811865475 ≈ 0.8408964152537145 b 1 = 1 × 0.7071067811865475 = 0.7071067811865475 ≈ 0.8408964152537145
t 1 = 0.25 − 1 × ( 1 − 0.8535533905932737 ) 2 = 0.25 − ( 0.1464466094067263 ) 2 t 1 = 0.25 − 1 × ( 1 − 0.8535533905932737 ) 2 = 0.25 − ( 0.1464466094067263 ) 2
t 1 = 0.25 − 0.0214466094067263 = 0.2285533905932737 t 1 = 0.25 − 0.0214466094067263 = 0.2285533905932737
p 1 = 2 × 1 = 2 p 1 = 2 × 1 = 2
π 1 ≈ ( 0.8535533905932737 + 0.8408964152537145 ) 2 4 × 0.2285533905932737 = 2.869662959118799 0.9142135623730948 ≈ 3.140579250522091 π 1 ≈ 4 × 0.2285533905932737 ( 0.8535533905932737 + 0.8408964152537145 ) 2 = 0.9142135623730948 2.869662959118799 ≈ 3.140579250522091
误差约 0.00101，有效数字约 3 位。
第 2 次迭代：
a 2 = 0.8535533905932737 + 0.8408964152537145 2 = 0.8472249029234941 a 2 = 2 0.8535533905932737 + 0.8408964152537145 = 0.8472249029234941
b 2 = 0.8535533905932737 × 0.8408964152537145 ≈ 0.8472012667468915 b 2 = 0.8535533905932737 × 0.8408964152537145 ≈ 0.8472012667468915
t 2 = 0.2285533905932737 − 2 × ( 0.8535533905932737 − 0.8472249029234941 ) 2 t 2 = 0.2285533905932737 − 2 × ( 0.8535533905932737 − 0.8472249029234941 ) 2
( a 1 − a 2 ) = 0.8535533905932737 − 0.8472249029234941 = 0.0063284876697796 ( a 1 − a 2 ) = 0.8535533905932737 − 0.8472249029234941 = 0.0063284876697796
( a 1 − a 2 ) 2 = 4.004975 × 10 − 5 ( a 1 − a 2 ) 2 = 4.004975 × 1 0 − 5
t 2 = 0.2285533905932737 − 2 × 4.004975 × 10 − 5 = 0.2284732910854960 t 2 = 0.2285533905932737 − 2 × 4.004975 × 1 0 − 5 = 0.2284732910854960
p 2 = 4 p 2 = 4
π 2 ≈ ( 0.8472249029234941 + 0.8472012667468915 ) 2 4 × 0.2284732910854960 = 2.870240247937589 0.9138931643419841 ≈ 3.141592646562934 π 2 ≈ 4 × 0.2284732910854960 ( 0.8472249029234941 + 0.8472012667468915 ) 2 = 0.9138931643419841 2.870240247937589 ≈ 3.141592646562934
误差约 7.1 × 10 − 9 7.1 × 1 0 − 9
第 3 次迭代：
a 3 = 0.8472249029234941 + 0.8472012667468915 2 = 0.8472130848351928 a 3 = 2 0.8472249029234941 + 0.8472012667468915 = 0.8472130848351928
b 3 = 0.8472249029234941 × 0.8472012667468915 ≈ 0.8472130847527654 b 3 = 0.8472249029234941 × 0.8472012667468915 ≈ 0.8472130847527654
t 3 = 0.2284732910854960 − 4 × ( 0.8472249029234941 − 0.8472130848351928 ) 2 t 3 = 0.2284732910854960 − 4 × ( 0.8472249029234941 − 0.8472130848351928 ) 2
( a 2 − a 3 ) = 0.8472249029234941 − 0.8472130848351928 = 1.18180883013 × 10 − 5 ( a 2 − a 3 ) = 0.8472249029234941 − 0.8472130848351928 = 1.18180883013 × 1 0 − 5
( a 2 − a 3 ) 2 = 1.39667 × 10 − 10 ( a 2 − a 3 ) 2 = 1.39667 × 1 0 − 10
t 3 = 0.2284732910854960 − 4 × 1.39667 × 10 − 10 = 0.2284732909488282 t 3 = 0.2284732910854960 − 4 × 1.39667 × 1 0 − 10 = 0.2284732909488282
p 3 = 8 p 3 = 8
π 3 ≈ ( 0.8472130848351928 + 0.8472130847527654 ) 2 4 × 0.2284732909488282 = 2.870241217692305 0.9138931637953128 ≈ 3.141592653589794 π 3 ≈ 4 × 0.2284732909488282 ( 0.8472130848351928 + 0.8472130847527654 ) 2 = 0.9138931637953128 2.870241217692305 ≈ 3.141592653589794
误差约 1.3 × 10 − 15 1.3 × 1 0 − 15
迭代次数 n n
π n π n 误差
有效数字位数
0
2.914213562373095
2.27 × 10 − 1 2.27 × 1 0 − 1 0-1
1
3.140579250522091
1.01 × 10 − 3 1.01 × 1 0 − 3 3
2
3.141592646562934
7.10 × 10 − 9 7.10 × 1 0 − 9 8-9
3
3.141592653589794
1.30 × 10 − 15 1.30 × 1 0 − 15 15-16
4
≈ ≈ < 10 − 30 < 1 0 − 30 30+
5
≈ ≈ < 10 − 60 < 1 0 − 60 60+
10
≈ ≈ < 10 − 2000 < 1 0 − 2000 2000+
这便是二次收敛 的威力：每步有效数字约翻倍。
为什么高斯-勒让德算法收敛如此之快？原因在于：
a n a n b n b n c n = a n − b n c n = a n − b n 平方速度 衰减具体来说，c n + 1 ≈ c n 2 8 a n c n + 1 ≈ 8 a n c n 2
我们可以从原理上证明：
第 n n ϵ n ≈ π n − π π ϵ n ≈ π π n − π ϵ n + 1 ≈ ϵ n 2 ϵ n + 1 ≈ ϵ n 2
设第 1 步后误差约 10 − 3 1 0 − 3
第 2 步：10 − 6 1 0 − 6 7.1 × 10 − 9 7.1 × 1 0 − 9
第 3 步：10 − 12 1 0 − 12
第 4 步：10 − 24 1 0 − 24
第 5 步：10 − 48 1 0 − 48
第 10 步：10 − 1536 1 0 − 1536
第 20 步：10 − 1572864 ≈ 10 − 1.5 × 10 6 1 0 − 1572864 ≈ 1 0 − 1.5 × 1 0 6
第 25 步：10 − 50331648 ≈ 10 − 5 × 10 7 1 0 − 50331648 ≈ 1 0 − 5 × 1 0 7
让我们看看高斯-勒让德算法与传统方法的收敛速度对比：
方法
收敛阶
每步新增位数
达到 1 亿位所需步数
阿基米德多边形法
线性
0.5 位
> 10 8 > 1 0 8
Machin 类反正切公式
线性
1-2 位
∼ 10 8 ∼ 1 0 8
高斯-勒让德算法 二次 翻倍 约 27 步
Chudnovsky 算法
超线性
14-16 位
约 700 万项
Borwein 四次算法
四次
翻四倍
约 14 步
注意：虽然 Chudnovsky 算法每步新增约 14 位，但它需要计算大量级数项 （每项需要高阶乘和幂运算）。相比之下，高斯-勒让德算法每步仅需计算一次平方根、一次加法和一次乘法，计算每步的代价远低于 Chudnovsky 。这正是高斯-勒让德算法在 1980-1990 年代占据 π π
from decimal import Decimal, getcontext
import math
def gauss_legendre_pi(precision_digits):
"""
使用高斯-勒让德算法计算 pi，精度为 precision_digits 位十进制小数。
参数:
precision_digits: 目标精度（十进制位数）
返回:
pi 的 Decimal 近似值
"""
# 设置 Decimal 精度（需要额外一些位数保证舍入正确）
getcontext().prec = precision_digits + 10
# 初始值
one = Decimal(1)
two = Decimal(2)
four = Decimal(4)
a = one # a0 = 1
b = one / Decimal(2).sqrt() # b0 = 1 / sqrt(2)
t = one / four # t0 = 1/4
p = one # p0 = 1
# 迭代次数估算：log2(precision_digits) + 1
# 因为每次迭代有效数字翻倍
iterations = precision_digits.bit_length() + 1
for _ in range(iterations):
a_next = (a + b) / two
b_next = (a * b).sqrt()
t = t - p * (a - a_next) ** 2
a, b = a_next, b_next
p = p * 2
# 计算最终 pi 近似
pi_approx = (a + b) ** 2 / (four * t)
return pi_approx
# 示例：计算 100 位 pi
pi_100 = gauss_legendre_pi(100)
print(f"Pi (100 位): {pi_100}")
# 与真实值比较
true_pi = "3.14159265358979323846264338327950288419716939937510" \
"58209749445923078164062862089986280348253421170679"
print(f"真实值: {true_pi}")
print(f"匹配: {str(pi_100)[:102] == true_pi}")
在 C++ 中使用 GMP（GNU Multiple Precision Arithmetic Library）或 MPFR 库实现：
#include <gmp.h>
#include <mpfr.h>
#include <cmath>
void gauss_legendre_pi(mpfr_t result, mpfr_prec_t prec) {
mpfr_t a, b, t, p, a_next, b_next, diff, temp;
// 初始化并分配内存
mpfr_inits2(prec, a, b, t, p, a_next, b_next, diff, temp, NULL);
// a0 = 1
mpfr_set_ui(a, 1, MPFR_RNDN);
// b0 = 1 / sqrt(2)
mpfr_sqrt_ui(b, 2, MPFR_RNDN);
mpfr_ui_div(b, 1, b, MPFR_RNDN);
// t0 = 0.25
mpfr_set_d(t, 0.25, MPFR_RNDN);
// p0 = 1
mpfr_set_ui(p, 1, MPFR_RNDN);
// 迭代次数: log2(prec) + 1
int iter = (int)log2(prec) + 1;
for (int i = 0; i < iter; i++) {
// a_next = (a + b) / 2
mpfr_add(a_next, a, b, MPFR_RNDN);
mpfr_div_ui(a_next, a_next, 2, MPFR_RNDN);
// b_next = sqrt(a * b)
mpfr_mul(b_next, a, b, MPFR_RNDN);
mpfr_sqrt(b_next, b_next, MPFR_RNDN);
// diff = a - a_next
mpfr_sub(diff, a, a_next, MPFR_RNDN);
// t = t - p * diff^2
mpfr_mul(diff, diff, diff, MPFR_RNDN);
mpfr_mul(diff, diff, p, MPFR_RNDN);
mpfr_sub(t, t, diff, MPFR_RNDN);
// a = a_next, b = b_next
mpfr_set(a, a_next, MPFR_RNDN);
mpfr_set(b, b_next, MPFR_RNDN);
// p = 2 * p
mpfr_mul_ui(p, p, 2, MPFR_RNDN);
}
// pi = (a + b)^2 / (4 * t)
mpfr_add(temp, a, b, MPFR_RNDN);
mpfr_mul(temp, temp, temp, MPFR_RNDN);
mpfr_mul_ui(t, t, 4, MPFR_RNDN);
mpfr_div(result, temp, t, MPFR_RNDN);
mpfr_clears(a, b, t, p, a_next, b_next, diff, temp, NULL);
}
操作
每步计算量
复杂度（对 n n
加法 a n + b n a n + b n
O ( n ) O ( n ) O ( n ) O ( n )
乘法 a n × b n a n × b n
O ( n log n ) O ( n log n ) O ( n log n ) O ( n log n )
平方根 x x
O ( n 2 ) O ( n 2 ) O ( n log n ) O ( n log n ) O ( n log n ) O ( n log n )
总计 O ( n log 2 n ) O ( n log 2 n )
利用 FFT（快速傅里叶变换）加速的大整数乘法，高斯-勒让德算法的整体复杂度为 O ( n log 2 n ) O ( n log 2 n ) π π
算法
收敛速度
每步复杂度
总复杂度
典型应用
高斯-勒让德 二次
O ( n log n ) O ( n log n ) O ( n log 2 n ) O ( n log 2 n ) 1980-1990s 记录保持者
Borwein 四次 四次
O ( n log n ) O ( n log n ) O ( n log 2 n ) O ( n log 2 n ) 类似性能，收敛更快
Chudnovsky 14-16 位/项
O ( n log n ) O ( n log n ) O ( n 2 ) O ( n 2 ) 当前百亿位记录
Machin 类 线性
O ( n ) O ( n ) O ( n 2 ) O ( n 2 ) 低精度快速计算
BBP 公式 线性
O ( n log n ) O ( n log n ) O ( n 2 ) O ( n 2 ) 提取任意十六进制位
高斯-勒让德算法在 π π
年份
记录位数
计算者
算法
1973
100 万
Guilloud & Bouyer
Machin 类公式
1976 200 万 Salamin 高斯-勒让德
1982
200 万
Tamura & Kanada
Machin 类公式
1983 1670 万 Kanada 等 高斯-勒让德
1989 5.3 亿 Chudnovsky 兄弟 Chudnovsky 算法
1999 2061 亿 Kanada 等 高斯-勒让德 (前 46 步)
2002
1.24 万亿
Kanada 等
双算法验证
2022
100 万亿
Google (Emma Haruka Iwao)
Chudnovsky
高斯-勒让德算法每步需要计算一个平方根——这是最耗时的操作之一。可以使用牛顿迭代法 加速：
对于计算 S S
x k + 1 = 1 2 ( x k + S x k ) x k + 1 = 2 1 ( x k + x k S )
由于高斯-勒让德的 a n ≈ b n a n ≈ b n
由于高精度计算极易出错，现代 π π 两种不同的算法同时计算 ，在最后比较结果以确保正确性。
常见的配对方案：
主算法
验证算法
优势
高斯-勒让德（二次）
Borwein 四次算法
共享 AGM 中间结果
Chudnovsky 级数
Machin 类公式
独立实现，验证可靠性
Chudnovsky 级数
AGM 类算法
跨家族验证
Kanada 的团队在 1999 年正是使用高斯-勒让德算法 和Borwein 四次算法 进行交叉验证，计算出了 2061 亿位 π π
高精度 π π
级别 1: 算法级并行
├── 主计算进程 (高斯-勒让德)
└── 验证进程 (Borwein 四次)
级别 2: 操作级并行
├── 大整数乘法 (FFT 并行)
├── 平方根计算 (牛顿迭代并行)
└── 加法和乘常数运算
级别 3: 数据级并行
├── 使用 SIMD 指令集加速 FFT
├── 多线程大整数运算
└── 分布式集群计算
高斯-勒让德算法的实际应用远超 π π
椭圆积分计算 ：AGM 是计算椭圆积分的最有效方法，广泛应用于天体力学、量子力学和工程领域对数和反三角函数的快速计算 ：利用 AGM 可以高效计算 ln x ln x arctan x arctan x 数值分析教学 ：是二次收敛算法的最佳教学案例之一高精度计算基准测试 ：新超级计算机常通过计算 π π
尽管高斯-勒让德算法有极快的收敛速度，它也存在一些局限：
局限
说明
数值稳定性 需要高精度算术，t n t n
内存需求 需要同时维护多个高精度数值
平方根瓶颈 每步的平方根运算在高精度下仍很昂贵
不可分离 不能像 BBP 公式那样直接提取 π π
舍入误差累积 多次舍入可能导致最终结果偏差
Borwein 兄弟（Jonathan 和 Peter Borwein）在 1984 年提出了四次收敛 的变体：
y 0 = 1 2 − 1 2 , a 0 = 6 − 4 2 y 0 = 2 1 − 2 1 , a 0 = 6 − 4 2
y n + 1 = 1 − ( 1 − y n 4 ) 1 / 4 1 + ( 1 − y n 4 ) 1 / 4 y n + 1 = 1 + ( 1 − y n 4 ) 1/4 1 − ( 1 − y n 4 ) 1/4
a n + 1 = a n ( 1 + y n + 1 ) 4 − 2 2 n + 3 y n + 1 ( 1 + y n + 1 + y n + 1 2 ) a n + 1 = a n ( 1 + y n + 1 ) 4 − 2 2 n + 3 y n + 1 ( 1 + y n + 1 + y n + 1 2 )
π n ≈ 1 a n （四次收敛，每步位数翻四倍） π n ≈ a n 1 （四次收敛，每步位数翻四倍）
Borwein 四次算法仅需约 17 步即可达到 100 亿位精度。
高斯-勒让德算法是数学计算史上的一座里程碑：
核心要点
说明
核心思想 利用算术-几何平均（AGM）迭代逼近 π π
收敛速度 二次收敛，每步有效数字翻倍
所需步骤 仅需约 log 2 ( n ) log 2 ( n ) n n
复杂度 O ( n log 2 n ) O ( n log 2 n )
历史意义 将 π π
现实应用 椭圆积分、初等函数计算、超算基准测试
一句话总结 ：高斯-勒让德算法通过算术-几何平均的二次收敛特性，以极少的迭代次数实现超高精度 π π
Brent, R. P. (1976). "Fast Multiple-Precision Evaluation of Elementary Functions". Journal of the ACM , 23(2), 242-251.
Salamin, E. (1976). "Computation of \pi Using Arithmetic-Geometric Mean". Mathematics of Computation , 30(135), 565-570.
Borwein, J. M. & Borwein, P. B. (1987). Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity . Wiley.
Bailey, D. H. (1988). "The Computation of \pi to 29,360,000 Decimal Digits Using Borweins' Quadratically Convergent Algorithm". Mathematics of Computation , 50(181), 283-296.
Kanada, Y. (2003). "Vectorization of the Gauss-Legendre Algorithm for Computing \pi". ACM Transactions on Mathematical Software , 29(3), 234-254.
Arndt, J. & Haenel, C. (2001). Pi - Unleashed . Springer-Verlag.