マチンの公式でどれくらいの精度で円周率を計算できるか？

次のようなグレゴリー級数で、$\pi/4$を計算することができます。

$\displaystyle \arctan(x) = \sum_{n=0}^{\infty}\frac{(-1)^n}{2n+1}x^{2n+1}=x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\frac{x^9}{9}\cdots$

そこで、この計算で100項までループした場合どれくらいの精度で計算できるか試してみました。

import decimal
d=lambda num:decimal.Decimal(str(num))
decimal.getcontext().prec = 10000

def arc(num,n):
    decimal.getcontext().prec = 10000
    total=0
    for i in range(n):
        x=2*i+1
        total+=(-1)**d(i)*(d(num)**d(x)/d(x))
    return total

pi_1=arc(d(1),100)*4
print(f'{pi_1:.20f}')

精度を出すために、decimal関数を使うことにしましたが、いちいちDecimalを打つのは面倒なので、ラムダ式を使いdという関数を作ってみました。また、arc関数では、arctanを求めたい値（num）と、第何項（n）まで求めるかを引数として渡すようにしました。答えは

3.13159290355855276431

驚くことに100項まで計算して、小数点第2位で誤差が発生してしまいます。大げさな割にはえらく情けない結果になってしまいました。そこで、少し工夫します。

$\displaystyle \tan(\alpha+\beta)=\frac{\tan \alpha+\tan \beta}{1-\tan \alpha\tan \beta}$

$\tan \alpha=a,\tan \beta=b$とすると

$\displaystyle \tan(\alpha+\beta)=\frac{a+b}{1-ab} $ から

$\displaystyle \alpha+\beta=\arctan\left(\frac{a+b}{1-ab}\right)$となります。

ここで$\displaystyle \frac{a+b}{1-ab}=1$となるa,bを探すと

$\displaystyle a=\frac{1}{2}; b=\frac{1}{3}$となります。

$\arctan 1=\displaystyle\arctan \frac{1}{2}+\arctan \frac{1}{3}$

それでは、Pythonで100項まで計算してみます。

p=100
machin=(arc(d(1)/d(2),p)+arc(d(1)/d(3),p))*4
print(machin)

小数点以下65桁まで求めると次の通りです。ちなみに次の行は正しい値です。
3.14159265358979323846264338327950288419716939937510582097494458734
3.14159265358979323846264338327950288419716939937510582097494459230
ということで小数点第61桁まで合っていることになります。1/2と1/3とarctanを2回計算するので手間は2倍になりますが、小数点第1位から61位まで一気に上がります。

このほか、arctanを使った式がたくさん見つかっています。特に有名なのがマチンの定理です。

$\displaystyle\frac{\pi}{4}=\displaystyle 4\arctan\frac{1}{5} – \arctan\frac{1}{239}=\displaystyle4\sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \left( \frac{1}{5}\right)^{2n+1}-\sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1}\left( \frac{1}{239}\right)^{2n+1} $

マチンの定理によると第100項まで求めて141桁まで行けてしまいます。さらに、木の発展型も沢山あります。

これを辞書形式にまとめ、何桁まで合っているか計算しました。結果は次の通りです。なお。桁数は次の記事のcompare関数を使いました。

円周率10万桁

def machin(name,parm,prec):
    total=0
    for i in parm:

        a=d(i[0])
        b=d(i[1])
        try:
            c=d(i[2])
        except IndexError:
            c=d(1)
        total+=a*arc(c/b,prec)
    return total*4
    
a_d={'arctan1':((1,1),),
     'arctan2':((1,2),(1,3)),
     'machin':((4,5),(-1,239)),
     'euler1':((5,7),(2,79,3)),
     'euler2':((4,5),(-1,70),(1,99)),
     'jacob':((2,2),(-1,7)),
     'hutton1':((3,4),(1,99,5)),
     'hutton2':((2,3),(1,7)),
     'gaus':((12,18),(8,57),(-5,239)),
     'stormer1':((6,8),(2,57),(1,239)),
     'stormer2':((44,57),(7,239),(-12,682),(24,12943)),
     'klingenstierna':((8,10),(-1,239),(-4,515)),
     'simson':((8,10),(-4,515),(-1,239)),
     'takano':((12,49),(32,57),(-5,239),(12,110443))}
for key,value in a_d.items():
    x=machin(key,value,100)
    print(f'{key:16}{x:.20f}  精度：{compare(pi_t,str(x))-2:3d}')

arctan1         3.13159290355855276431  精度：  1
arctan2         3.14159265358979323846  精度： 61
machin          3.14159265358979323846  精度：141
euler1          3.14159265358979323846  精度：170
euler2          3.14159265358979323846  精度：141
jacob           3.14159265358979323846  精度： 61
hutton1         3.14159265358979323846  精度：121
hutton2         3.14159265358979323846  精度： 97
gaus            3.14159265358979323846  精度：252
stormer1        3.14159265358979323846  精度：182
stormer2        3.14159265358979323846  精度：352
klingenstierna  3.14159265358979323846  精度：201
simson          3.14159265358979323846  精度：201
takano          3.14159265358979323846  精度：340

結果は上記の通りです。神奈川県の学校の教諭である高野喜久雄さんがさんが導き出した方法がガウスやオイラーを上回った結果を出しています。この公式を使って日立のスーパーコンピュータを円周率を1兆2411億桁まで計算したそうです。