# Copyright 2018 The go-python Authors. All rights reserved.

# Use of this source code is governed by a BSD-style

# license that can be found in the LICENSE file.

"""

Python3 program to calculate Pi using python long integers, binary

splitting and the Chudnovsky algorithm

See: http://www.craig-wood.com/nick/articles/ FIXME for explanation

Nick Craig-Wood <nick@craig-wood.com>

"""

import math

from time import time

def sqrt(n, one):

"""

Return the square root of n as a fixed point number with the one

passed in. It uses a second order Newton-Raphson convgence. This

doubles the number of significant figures on each iteration.

"""

# Use floating point arithmetic to make an initial guess

floating_point_precision = 10**16

n_float = float((n * floating_point_precision) // one) / floating_point_precision

x = (int(floating_point_precision * math.sqrt(n_float)) * one) // floating_point_precision

n_one = n * one

while 1:

x_old = x

x = (x + n_one // x) // 2

if x == x_old:

break

return x

def pi_chudnovsky_bs(digits):

"""

Compute int(pi * 10**digits)

This is done using Chudnovsky's series with binary splitting

"""

C = 640320

C3_OVER_24 = C**3 // 24

def bs(a, b):

"""

Computes the terms for binary splitting the Chudnovsky infinite series

a(a) = +/- (13591409 + 545140134*a)

p(a) = (6*a-5)*(2*a-1)*(6*a-1)

b(a) = 1

q(a) = a*a*a*C3_OVER_24

returns P(a,b), Q(a,b) and T(a,b)

"""

if b - a == 1:

# Directly compute P(a,a+1), Q(a,a+1) and T(a,a+1)

if a == 0:

Pab = Qab = 1

else:

Pab = (6*a-5)*(2*a-1)*(6*a-1)

Qab = a*a*a*C3_OVER_24

Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a)

if a & 1:

Tab = -Tab

else:

# Recursively compute P(a,b), Q(a,b) and T(a,b)

# m is the midpoint of a and b

m = (a + b) // 2

# Recursively calculate P(a,m), Q(a,m) and T(a,m)

Pam, Qam, Tam = bs(a, m)

# Recursively calculate P(m,b), Q(m,b) and T(m,b)

Pmb, Qmb, Tmb = bs(m, b)

# Now combine

Pab = Pam * Pmb

Qab = Qam * Qmb

Tab = Qmb * Tam + Pam * Tmb

return Pab, Qab, Tab

# how many terms to compute

DIGITS_PER_TERM = math.log10(C3_OVER_24/6/2/6)

N = int(digits/DIGITS_PER_TERM + 1)

# Calclate P(0,N) and Q(0,N)

P, Q, T = bs(0, N)

one = 10**digits

sqrtC = sqrt(10005*one, one)

return (Q*426880*sqrtC) // T

# The last 5 digits or pi for various numbers of digits

check_digits = (

(100, 70679),

(1000, 1989),

(10000, 75678),

(100000, 24646),

(1000000, 58151),

(10000000, 55897),

)

if __name__ == "__main__":

digits = 100

pi = pi_chudnovsky_bs(digits)

print(str(pi))

for digits, check in check_digits:

start =time()

pi = pi_chudnovsky_bs(digits)

print("chudnovsky_gmpy_bs: digits",digits,"time",time()-start)

last_five_digits = pi % 100000

if check == last_five_digits:

print("Last 5 digits %05d OK" % last_five_digits)

else:

print("Last 5 digits %05d wrong should be %05d" % (last_five_digits, check))