pi, . ""; mpmath , "" pi e. , Decimal , .
''' The Salamin / Brent / Gauss Arithmetic-Geometric Mean pi formula.
Let A[0] = 1, B[0] = 1/Sqrt(2)
Then iterate from 1 to 'n'.
A[n] = (A[n-1] + B[n-1])/2
B[n] = Sqrt(A[n-1]*B[n-1])
C[n] = (A[n-1]-B[n-1])/2
PI[n] = 4A[n+1]^2 / (1-(Sum (for j=1 to n; 2^(j+1))*C[j]^2))
Written by PM 2Ring 2008.10.19
Converted to use Decimal 2014.10.21
'''
import sys
from decimal import Decimal, getcontext, ROUND_DOWN
def AGM_pi(m):
a, b = Decimal(1), Decimal(2).sqrt() / 2
s, k = Decimal(0), Decimal(4)
for i in xrange(m):
c = (a - b) / 2
a, b = (a + b) / 2, (a * b).sqrt()
s += k * c * c
if False:
pi = 4 * a * a / (1 - s)
print "%2d:\n%s\n" % (i, pi)
k *= 2
return 4 * a * a / (1 - s)
def main():
prec = len(sys.argv) > 1 and int(sys.argv[1]) or 50
delta = 3
prec += 1 + delta
ctx = getcontext()
ctx.prec = prec
i, m = prec, 0
while i>0:
i >>= 1; m += 1
pi = AGM_pi(m)
ctx.rounding = ROUND_DOWN
ctx.prec -= delta
pi = +pi
print "pi ~=\n%s" % pi
if __name__ == '__main__':
main()
source
share