
Multiply one side of a hexagon by the radius (of its circumcircle), then multiply this by three, to yield the area of a dodecagon; if we cut a hexagon into a dodecagon, multiply its side by its radius, then again multiply by six, we get the area of a 24-gon; the finer we cut, the smaller the loss with respect to the area of circle, thus with further cut after cut, the area of the resulting polygon will coincide and become one with the circle; there will be no loss
Alright, so you start with a hexagon. This is approximately a circle, although only crudely. Divide the circumference by the radius, you approximate pi! Then, you subdivide the number of sides by two, so a hexagon becomes a dodecagon. Do the same thing, again and again!
Easy enough. Then, you really only need a robust function to calculate square roots and you're golden. I'll cheat and use ChatGPT, but let's solve this from scratch in Python (I know, ungodly slow for this kind of thing, but it's easier to read):
Code:
from decimal import Decimal, getcontext
def sqrt(x):
x = Decimal(x)
y, prev_y = x / 2, 0
while y != prev_y:
prev_y, y = y, (y + x / y) / 2
return y
def liu_hui_pi(n):
m, a = Decimal(6), Decimal(1)
for _ in range(n):
a = sqrt(Decimal(2) - sqrt(Decimal(4) - a**2))
m *= 2
return (m * a) / 2
decimal.getcontext().prec = 128
for n in range(2, 12):
print(f"{n} = {liu_hui_pi(n)}")
Accurate digits at iteration 25: 3.14159265358979311
Pretty sweet.
Only half of the digits will be accurate with respect to the Decimal precision, because of the square operation.
ADDENDUM
This is why I'm not a mathematician! They have calculated 202 trillion digits of pi using this method:
Chudnovsky algorithm
Code:
import decimal
def binary_split(a, b):
if b == a + 1:
Pab = -(6*a - 5)*(2*a - 1)*(6*a - 1)
Qab = 10939058860032000 * a**3
Rab = Pab * (545140134*a + 13591409)
else:
m = (a + b) // 2
Pam, Qam, Ram = binary_split(a, m)
Pmb, Qmb, Rmb = binary_split(m, b)
Pab = Pam * Pmb
Qab = Qam * Qmb
Rab = Qmb * Ram + Pam * Rmb
return Pab, Qab, Rab
def chudnovsky(n):
P1n, Q1n, R1n = binary_split(1, n)
return (426880 * decimal.Decimal(10005).sqrt() * Q1n) / (13591409*Q1n + R1n)
decimal.getcontext().prec = 128
for n in range(2, 12):
print(f"{n} = {chudnovsky(n)}")
Bookmarks