gistfile1.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

# マチンの公式
# π /4 = 4arctan(1/5) - arctan(1/239)
# π = (16/5 - 16/(3*5^3) + 16/(5*5^5) - ...) - (4/239 - 4/(3*239^3) + ...)
#   = (16/5 - 4/239) - (16/(3*5^3) - 4/(3*239^3)) + ...
# w[0] = 16/5,  w[1] = w0/(3*5^3),   w[2] = w1/(5*5^5), ...
# v[0] = 4/239, v[1] = v0/(3*239^3), v[2] = v1/(5*239^5), ...
# m[n] = ± (w[n-1]/(2n-1)(5^2) - v[n-1]/(2n-1)(239^2))
#      = ± (w[n-1]/5^2 - v[n-1]/239^2)/(2n-1) (n % 2 == 1 の時に+)

# ガウス＝ルジャンドルのアルゴリズム
# Library dependence:
#   gmpy
# 初期値
#   a[0] = 1, b[0] = 1 / √2, t[0] = 1 / 4, p[0] = 1
# 漸化式
#   a[n+1] = (a[n] + b[n]) / 2
#   b[n+1] = √(a[n] * b[n])
#   t[n+1] = t[n] - p[n] * (a[n] - a[n+1])^2
#   p[p+1] = 2 * p[n]

import math

# デフォルト計算桁数
_L = 1000

# マチンの公式
# Machin-like formula
def machin(l = _L):
    d = (l / 4) + 1
    n = (l / (2 * math.log(5, 10))) + 1

    w = 16 * 5 * (10 ** l)
    v = 4 * 239 * (10 ** l)
    divw = 5 ** 2
    divv = 239 ** 2

    p = 0
    for i in xrange(1, int(n)):
        w = w / divw
        v = v / divv
        m = (w - v) / (2 * i - 1)
        p = p + (m if i % 2 else -m)
    return p


# ガウス＝ルジャンドルのアルゴリズム
# Gauss–Legendre algorithm
# Dependent(依存): gmpy
def gauss_legendre(l = _L):
    import gmpy
    bits = int(l / math.log(2, 10)) + 1
    # 一度のループごとに 2^n 桁以上に収束していく
    n = int((math.log(l, 10) / math.log(2, 10)) + 1)
    a = gmpy.mpf(1, bits)
    b = gmpy.fsqrt(gmpy.mpf(2, bits)) / 2
    t = gmpy.mpf(1.0 / 4, bits)
    p = gmpy.mpz(1)

    # 何回計算すれば目的の桁に達するのかを得る
    for _ in xrange(n):
        a_1 = a  # a[n-1]
        a = (a + b) / 2
        b = gmpy.fsqrt(a_1 * b)
        t = t - p * ((a_1 - a) ** 2)
        p = p * 2
    
    return ((a + b) ** 2) / (4 * t)


def main():
    print '',machin()
    print gauss_legendre()


if __name__ == '__main__':
    main()