|
| 1 | +""" |
| 2 | +Python3 program to calculate Pi using python long integers, binary |
| 3 | +splitting and the Chudnovsky algorithm |
| 4 | +
|
| 5 | +See: http://www.craig-wood.com/nick/articles/ FIXME for explanation |
| 6 | +
|
| 7 | +Nick Craig-Wood <nick@craig-wood.com> |
| 8 | +""" |
| 9 | + |
| 10 | +import math |
| 11 | +from time import time |
| 12 | + |
| 13 | +def sqrt(n, one): |
| 14 | + """ |
| 15 | + Return the square root of n as a fixed point number with the one |
| 16 | + passed in. It uses a second order Newton-Raphson convgence. This |
| 17 | + doubles the number of significant figures on each iteration. |
| 18 | + """ |
| 19 | + # Use floating point arithmetic to make an initial guess |
| 20 | + floating_point_precision = 10**16 |
| 21 | + n_float = float((n * floating_point_precision) // one) / floating_point_precision |
| 22 | + x = (int(floating_point_precision * math.sqrt(n_float)) * one) // floating_point_precision |
| 23 | + n_one = n * one |
| 24 | + while 1: |
| 25 | + x_old = x |
| 26 | + x = (x + n_one // x) // 2 |
| 27 | + if x == x_old: |
| 28 | + break |
| 29 | + return x |
| 30 | + |
| 31 | + |
| 32 | +def pi_chudnovsky_bs(digits): |
| 33 | + """ |
| 34 | + Compute int(pi * 10**digits) |
| 35 | +
|
| 36 | + This is done using Chudnovsky's series with binary splitting |
| 37 | + """ |
| 38 | + C = 640320 |
| 39 | + C3_OVER_24 = C**3 // 24 |
| 40 | + def bs(a, b): |
| 41 | + """ |
| 42 | + Computes the terms for binary splitting the Chudnovsky infinite series |
| 43 | +
|
| 44 | + a(a) = +/- (13591409 + 545140134*a) |
| 45 | + p(a) = (6*a-5)*(2*a-1)*(6*a-1) |
| 46 | + b(a) = 1 |
| 47 | + q(a) = a*a*a*C3_OVER_24 |
| 48 | +
|
| 49 | + returns P(a,b), Q(a,b) and T(a,b) |
| 50 | + """ |
| 51 | + if b - a == 1: |
| 52 | + # Directly compute P(a,a+1), Q(a,a+1) and T(a,a+1) |
| 53 | + if a == 0: |
| 54 | + Pab = Qab = 1 |
| 55 | + else: |
| 56 | + Pab = (6*a-5)*(2*a-1)*(6*a-1) |
| 57 | + Qab = a*a*a*C3_OVER_24 |
| 58 | + Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a) |
| 59 | + if a & 1: |
| 60 | + Tab = -Tab |
| 61 | + else: |
| 62 | + # Recursively compute P(a,b), Q(a,b) and T(a,b) |
| 63 | + # m is the midpoint of a and b |
| 64 | + m = (a + b) // 2 |
| 65 | + # Recursively calculate P(a,m), Q(a,m) and T(a,m) |
| 66 | + Pam, Qam, Tam = bs(a, m) |
| 67 | + # Recursively calculate P(m,b), Q(m,b) and T(m,b) |
| 68 | + Pmb, Qmb, Tmb = bs(m, b) |
| 69 | + # Now combine |
| 70 | + Pab = Pam * Pmb |
| 71 | + Qab = Qam * Qmb |
| 72 | + Tab = Qmb * Tam + Pam * Tmb |
| 73 | + return Pab, Qab, Tab |
| 74 | + # how many terms to compute |
| 75 | + DIGITS_PER_TERM = math.log10(C3_OVER_24/6/2/6) |
| 76 | + N = int(digits/DIGITS_PER_TERM + 1) |
| 77 | + # Calclate P(0,N) and Q(0,N) |
| 78 | + P, Q, T = bs(0, N) |
| 79 | + one = 10**digits |
| 80 | + sqrtC = sqrt(10005*one, one) |
| 81 | + return (Q*426880*sqrtC) // T |
| 82 | + |
| 83 | +# The last 5 digits or pi for various numbers of digits |
| 84 | +check_digits = ( |
| 85 | + (100, 70679), |
| 86 | + (1000, 1989), |
| 87 | + (10000, 75678), |
| 88 | + (100000, 24646), |
| 89 | + (1000000, 58151), |
| 90 | + (10000000, 55897), |
| 91 | +) |
| 92 | + |
| 93 | +if __name__ == "__main__": |
| 94 | + digits = 100 |
| 95 | + pi = pi_chudnovsky_bs(digits) |
| 96 | + print(str(pi)) |
| 97 | + for digits, check in check_digits: |
| 98 | + start =time() |
| 99 | + pi = pi_chudnovsky_bs(digits) |
| 100 | + print("chudnovsky_gmpy_bs: digits",digits,"time",time()-start) |
| 101 | + last_five_digits = pi % 100000 |
| 102 | + if check == last_five_digits: |
| 103 | + print("Last 5 digits %05d OK" % last_five_digits) |
| 104 | + else: |
| 105 | + print("Last 5 digits %05d wrong should be %05d" % (last_five_digits, check)) |
0 commit comments