Skip to content

Instantly share code, notes, and snippets.

@PM2Ring
Created October 31, 2023 18:39
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save PM2Ring/032ffa274fec417c47cc31f014dcfa5c to your computer and use it in GitHub Desktop.
Save PM2Ring/032ffa274fec417c47cc31f014dcfa5c to your computer and use it in GitHub Desktop.
The Salamin / Brent / Gauss Arithmetic-Geometric Mean pi formula, using fixed point arithmetic
""" Gauss / Salamin / Brent true AGM pi
Binary fixed point
Written by PM 2Ring 2023.11.1
"""
from math import isqrt
def AGM_pi(loops, bits):
a, b = 1 << bits, isqrt(1 << 2*bits-1)
s = 1 << bits-2
for i in range(loops):
c = (a - b) >> 1
s -= c * c >> bits-i
a, b = (a + b) >> 1, isqrt(a * b)
return a * a // s
# Tested upto 18 loops
loops = 11
print(f"{loops=}")
bits = ((161139 << loops) - 165914) // 35553
print(f"{bits=}")
digits = bits * 30103 // 100000
print(f"{digits=}")
prec = bits + 8 + loops
print(f"{prec=}")
mypi = AGM_pi(loops, prec) * 10**digits >> prec
#print(mypi)
mypi = str(mypi)[1:]
print("pi ~=\n3.")
for i in range(0, len(mypi), 5):
ii = i + 5
j = 0 if ii % 25 else 1 if ii % 100 else 2 if ii % 1000 else 3
print(mypi[i:ii], end=j*"\n" or " ")
@PM2Ring
Copy link
Author

PM2Ring commented Oct 31, 2023

@PM2Ring
Copy link
Author

PM2Ring commented Oct 31, 2023

bits = (2**loops * pi - log(8*pi)) / log(2)
See (PDF) Jonathan Borwein, Pi and the AGM by Richard P. Brent.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment