The Gauss–Brent–Salamin formula for π

Gauss's theorem expresses K(k) by the AGM M(1, k′), the AGM evaluation of the second kind expresses E(k) / K(k) by the squared AGM corrections, and Legendre's relation ties the complete integrals at complementary moduli k and k′ = √(1 - k²) to π / 2. Combined, they express π by two AGM iterations alone:

• TauCeti.Real.pi_eq_four_mul_agm_mul_agm_div_one_sub_tsum: for x, y > 0 with x² + y² = 1, Salamin's formula π = 4 M(1, x) M(1, y) / (1 - ∑ 2ʲ (cⱼ(x)² + cⱼ(y)²));
• TauCeti.Real.pi_eq_four_mul_agm_sq_div_one_sub_tsum: the self-complementary case x = y = 1 / √2, the Gauss–Brent–Salamin formula π = 4 M(1, 1/√2)² / (1 - ∑_{j ≥ 1} 2^(j+1) cⱼ²).

Here cⱼ² = aⱼ² - gⱼ² for the AGM iteration a₀ = 1, g₀ = x, with j ≥ 1. Mathlib's NNReal.agmSequences 1 x n is the pair (gₙ₊₁, aₙ₊₁), so its index n corresponds to the classical index j = n + 1. Truncating the series and the AGM after finitely many steps gives the Brent–Salamin algorithm for computing π.

References

• E. Salamin, Computation of π using arithmetic-geometric mean, Math. Comp. 30 (1976), 565–570.
• R. P. Brent, Fast multiple-precision evaluation of elementary functions, J. ACM 23 (1976), 242–251.
• J. M. Borwein and P. B. Borwein, Pi and the AGM, Wiley, 1987, §2.4.

theorem TauCeti.Real.pi_eq_four_mul_agm_mul_agm_div_one_sub_tsum {x y : NNReal} (hx : 0 < x) (hy : 0 < y) (hxy : x ^ 2 + y ^ 2 = 1) : Real.pi = 4 * ↑(NNReal.agm 1 x) * ↑(NNReal.agm 1 y) / (1 - ∑' (n : ℕ), 2 ^ (n + 1) * (↑(NNReal.agmSequences 1 x n).2 ^ 2 - ↑(NNReal.agmSequences 1 x n).1 ^ 2 + (↑(NNReal.agmSequences 1 y n).2 ^ 2 - ↑(NNReal.agmSequences 1 y n).1 ^ 2)))

Salamin's formula for complementary moduli: if x, y > 0 and x² + y² = 1, then π = 4 M(1, x) M(1, y) / (1 - ∑ₙ 2ⁿ⁺¹ (cₙ₊₁(x)² + cₙ₊₁(y)²)), where M is the AGM and cₙ₊₁(x)² = aₙ₊₁² - gₙ₊₁² is the squared correction of the AGM iteration started at (1, x). Mathlib's agmSequences 1 x n is the pair (gₙ₊₁, aₙ₊₁).

theorem TauCeti.Real.pi_eq_four_mul_agm_sq_div_one_sub_tsum : Real.pi = 4 * ↑(NNReal.agm 1 (NNReal.sqrt 2)⁻¹) ^ 2 / (1 - ∑' (n : ℕ), 2 ^ (n + 2) * (↑(NNReal.agmSequences 1 (NNReal.sqrt 2)⁻¹ n).2 ^ 2 - ↑(NNReal.agmSequences 1 (NNReal.sqrt 2)⁻¹ n).1 ^ 2))

The Gauss–Brent–Salamin formula: with the AGM iteration a₀ = 1, g₀ = 1 / √2 and cⱼ² = aⱼ² - gⱼ², π = 4 M(1, 1/√2)² / (1 - ∑_{j ≥ 1} 2^(j+1) cⱼ²). Mathlib's agmSequences 1 (√2)⁻¹ n is the pair (gₙ₊₁, aₙ₊₁), so the series index is n = j - 1.