● 各種 atan() 展開式の演算量比較

　もう一度 atan() の展開式をながめます。
　atan(1/n) = 1 / (1*n) - 1 / (3*n³) + 1 / (5*n⁵) - 1 / (7*n⁷) + ...

　この右辺は 1 / (m * n^m) の集計と考えることができます。(m = 1, 3, 5, ... )
おおざっばに言って、ある n について m を大きくしていって、1 / (m * n^m) が限度を超えて小さくなれば、以後の集計処理は無駄でしょう。 atan の値を集計記憶している変数に収まりきらなくなるからです。 限度とは、5000 桁求めるのであれば 1e-5000 程度です。

　各種の atan() によるπ展開式について、演算量を見積もるプログラムを作りました。

　ソースプログラム

　実行結果：

n     2, m 16597, eps 1e-5000.41
n     5, m  7149, eps 1e-5000.79
n     8, m  5533, eps 1e-5000.54

n     2, m 16597, eps 1e-5000.41
n     3, m 10473, eps 1e-5000.91

n     2, m 16597, eps 1e-5000.41
n     7, m  5913, eps 1e-5000.84

n     3, m 10473, eps 1e-5000.91
n     7, m  5913, eps 1e-5000.84

n     4, m  8299, eps 1e-5000.41
n    20, m  3841, eps 1e-5000.84
n  1985, m  1517, eps 1e-5005.88

n     5, m  7149, eps 1e-5000.79
n    70, m  2709, eps 1e-5001.8
n    99, m  2505, eps 1e-5002.46

n     5, m  7149, eps 1e-5000.79
n   408, m  1915, eps 1e-5002.7
n  1393, m  1591, eps 1e-5005.23

n     5, m  7149, eps 1e-5000.79
n   239, m  2101, eps 1e-5000.34

n    18, m  3981, eps 1e-5000.84
n    57, m  2847, eps 1e-5002.43
n   239, m  2101, eps 1e-5000.34

 Name        m1     m2     m3            total
-----------------------------------------------
Dase       16597   7149   5533      0 :  29279
Euler      16597  10473      0      0 :  27070
Vega       16597   5913      0      0 :  22510
Clausen    10473   5913      0      0 :  16386
Gauss (1)   8299   3841   1517      0 :  13657
Rutherford  7149   2709   2505      0 :  12363
Vega        7149   1915   1591      0 :  10655
Machin      7149   2101      0      0 :   9250
Gauss (2)   3981   2847   2101      0 :   8929

　Name は展開式の発見者の名前、m1 〜 m3 は atan(1/n) それぞれ項について、m をどれほど大きくしたら限度を超えられそうかの見積り、total は m1 〜 m3 の合計です。

　n が大きくなると m はそれほど大きくしなくてもよくなりますが、n が 100 あたりを越えると、影響が小さくなります。 マチンの展開式は項数が少ない割に演算量が少ないため、よく使われるのでしょう。

　前記のプログラムは 1 / (m * n^m) < limit となる m を求めるのに何も工夫していません。 ニュートン・ラプソン法で m を求めるのなら、次のようにします。 pm(5, log(1e-300)) のように使います。

　#include <stdio.h>
　#include <math.h>
　#include <stdint.h>
　/*
　    f(x) = log(x) + a * x + b, f(x) = 0
　    f'(x) = 1/x + a
　    x(n+1) = x(n) - f(x(n))/f'(x(n))
　*/
　int pm(int n, double loglimit) {
　    int i;
　    double x, fx, dx, logn;
　    logn = log(n);
　    x = 33;　　　　　  /* 適当な初期値 */
　    for (i = 0; i < 20; i += 1) {
　        fx = log(x) + logn * x + loglimit;
　        dx = fx / (1 / x + logn);
　        if (fabs(dx / x) < 1e-17)
　            break;
　        x -= dx;
　    }
　    int m = (i < 20) ? lround(x) : INT32_MIN;
　    if ((m & 1) == 0)
　        m += 1;        /* 条件を満足する最小の奇数 */
　    return m;
　}