pi1000.py

import copy
import typing

def sign(x: typing.Union[int, float, 'Decimal1000']) -> int:
    """x の符号が正なら 1 、負なら -1 、そうでないなら 0 を返す"""
    if isinstance(x, Decimal1000):
        return x.sign()

    if x < 0: return -1
    if x > 0: return 1
    return 0

NUM_CHUNKS = 252
NUM_SIGNIFICANT_CHUNKS = 251
MOD_MASK = 10000

class Decimal1000:
    """整数 + 小数点以下 1000 桁"""

    # 各種演算の実装方針
    # _chunks を毎回コピーしていたらつらいので、二項算術演算子は実装せずミュータブルな構造にする

    def __init__(self, x: typing.Union[int, 'Decimal1000']) -> None:
        if isinstance(x, int):
            # _chunks[0] → 整数部
            # _chunks[1~250] → 小数部（4桁ずつ）
            # _chunks[251] → ×16 とかの操作があるので、精度に余裕を持たせる
            self._chunks = [0] * NUM_CHUNKS
            self._chunks[0] = abs(x)
            self._sign = sign(x)
        elif isinstance(x, Decimal1000):
            self.assign(x)
        else:
            raise TypeError()

    def assign(self, x: 'Decimal1000') -> None:
        """x を self にコピー"""
        if hasattr(self, '_chunks'):
            for i in range(NUM_CHUNKS):
                self._chunks[i] = x._chunks[i]
        else:
            self._chunks = copy.copy(x._chunks)
        self._sign = x._sign

    def __copy__(self) -> 'Decimal1000':
        return Decimal1000(self)

    __deepcopy__ = __copy__

    def sign(self) -> int:
        return self._sign

    def iszero(self) -> bool:
        return self._sign == 0

    def isnegative(self) -> bool:
        return self._sign < 0

    def negate(self) -> None:
        """符号を反転"""
        self._sign = -self._sign

    def setzero(self) -> None:
        """self = 0"""
        for chunk_index in range(NUM_CHUNKS):
            self._chunks[chunk_index] = 0
        self._sign = 0

    def add(self, x: typing.Union[int, 'Decimal1000']) -> None:
        """self += x"""
        if isinstance(x, int):
            if x != 0:
                # 整数部で足し算
                i = (self._sign * self._chunks[0]) + x
                self._chunks[0] = abs(i)
                self._sign = sign(i)
        elif isinstance(x, Decimal1000):
            if x.iszero():
                pass
            elif self.iszero():
                # self = 0 ならコピーしてくるだけ
                self.assign(x)
            elif self._sign == x._sign:
                # 符号が一致しているので足し合わせるだけ
                self._add(self._chunks, x._chunks)
            elif x.isnegative() and not self.isnegative():
                # self = self - |x|
                self._sub(self._chunks, x._chunks)
            else:
                # self = x - |self|
                self._sub(x._chunks, self._chunks)
        else:
            raise TypeError()

    def _add(self, x: typing.List[int], y: typing.List[int]) -> None:
        """self = x + y"""
        carry = 0

        # 小数部の計算
        for chunk_index in range(NUM_CHUNKS - 1, 0, -1):
            s = x[chunk_index] + y[chunk_index] + carry
            self._chunks[chunk_index] = s % MOD_MASK
            carry = s // MOD_MASK

        # 整数部の計算
        self._chunks[0] = x[0] + y[0] + carry

    def _sub(self, x: typing.List[int], y: typing.List[int]) -> None:
        """self = x - y"""

        for chunk_index in range(NUM_CHUNKS):
            if x[chunk_index] > y[chunk_index]:
                # x > y
                self._sub_core(x, y)
                self._sign = 1
                return
            if x[chunk_index] < y[chunk_index]:
                # x < y
                self._sub_core(y, x)
                self._sign = -1
                return

        # x = y なので self = 0
        self.setzero()

    def _sub_core(self, x: typing.List[int], y: typing.List[int]) -> None:
        """self = x - y (x > y)"""
        borrow = False

        # 小数部の計算
        for chunk_index in range(NUM_CHUNKS - 1, 0, -1):
            left = x[chunk_index]
            right = y[chunk_index]

            # 前回の桁借りを反映
            if borrow: left -= 1

            # 新たに桁借りが必要か？
            borrow = left < right
            if borrow: left += MOD_MASK

            self._chunks[chunk_index] = left - right

        # 整数部の計算
        left = x[0]
        right = y[0]
        if borrow: left -= 1
        assert left >= right
        self._chunks[0] = left - right

    def mul(self, x: int) -> None:
        """self *= x"""
        if not isinstance(x, int): raise TypeError()
        if self.iszero(): return

        if x == 0:
            self.setzero()
            return

        if x < 0:
            self.negate()
            x = -x

        carry = 0

        # すべてのチャンクに x を掛ける
        for chunk_index in range(NUM_CHUNKS - 1, 0, -1):
            p = self._chunks[chunk_index] * x + carry
            self._chunks[chunk_index] = p % MOD_MASK
            carry = p // MOD_MASK

        self._chunks[0] = self._chunks[0] * x + carry

    def div(self, x: int) -> None:
        """self /= x"""
        if not isinstance(x, int): raise TypeError()
        if x == 0: raise ZeroDivisionError()
        if self.iszero(): return

        if x < 0:
            self.negate()
            x = -x

        # x で割って、余りを次に持ち越す
        # 一番最後の桁は切り捨て
        r = 0
        for chunk_index in range(NUM_CHUNKS):
            left = self._chunks[chunk_index] + (r * MOD_MASK)
            q = left // x
            if chunk_index >= 1: assert q < MOD_MASK
            self._chunks[chunk_index] = q
            r = left % x

    def __str__(self) -> str:
        """小数点以下 1000 桁までを文字列として返す"""
        s = str(self._chunks[0]) + '.'
        for i in range(1, NUM_SIGNIFICANT_CHUNKS):
            s += '{:04d}'.format(self._chunks[i])
        return s

def highres_arctan_inv(y: int, n: int) -> Decimal1000:
    # P_0
    p = Decimal1000(y)
    p.div(y * y)
    # C_0
    c = 1
    # T_0
    t = Decimal1000(p)
    t.div(c)
    s = Decimal1000(t)

    for i in range(1, n + 1):
        # P_i = P_i-1 / (y * y)
        p.div(y * y)
        # C_i = C_i-1 + 2
        c += 2

        # T_i = (-1)^i * (P_i / C_i)
        t.assign(p)
        t.div(c)
        if i % 2 != 0:
            t.negate()

        s.add(t)

    return s

def main() -> None:
    # M1 = 16 * arctan(1 / 5)
    pi = highres_arctan_inv(5, 714)
    pi.mul(16)

    # M2 = -4 * arctan(1 / 239)
    m2 = highres_arctan_inv(239, 210)
    m2.mul(-4)

    # π = M1 + M2
    pi.add(m2)

    print(pi)

    # 正しい結果: http://www.geocities.jp/f9305710/PAI1000000.html
    ans = '3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989'
    print('チェック: ' + ('OK' if str(pi) == ans else 'NG'))

if __name__ == '__main__':
    main()