円周率.jp (http://xn--w6q13e505b.jp/method/binary.html)

Binary Splitting

Binary Splitting Method

あえて日本語で言うなら「二分割法」とでもいうかもしれないが、 現状日本語での言及も「Binary Splitting」や「BS 法」 と言及されるばかりである。 初出の論文である [FT05] によると各関数や定数についていえば Karatsuba が同様の概念を提案していたらしい?

$n$ 桁乗算の計算量を $M(n)$ とするとき、分割統治を行うことで 有理級数が $O(\log n\cdot M(n\log n))$ の計算量で計算できる手法である。 [FT06] 現在 $M(n)=O(n\log n)$ とされているので、$O(n(\log n)^3)$ になる。

簡易版

まずは広く使われている形式を紹介する。 計算しようとする対象が \[ \frac{P(0,N)}{Q(0,N)} = \sum_{n=0}^{N-1} \frac{a(n)}{b(n)} \frac{p(0)\cdots p(n)}{q(0)\cdots q(n)} \]

のような形になるよう、 適切な式 $P(l,u)$、$Q(l,u)$、$R(l,u)$ を設定する。 (具体例は下記) するとそれらの具体的計算については $l\le m\lt u$ となる $m$ を用いて以下のように再帰計算ができる。 $l+1=u$ となる場合は定義式($a(n)$ など)に $n=u$ を代入する。

\[ \begin{eqnarray} P(l, u) &=& P(l, m)Q(m, u) + R(l, m)P(m, u)\\ Q(l, u) &=& Q(l, m)Q(m, u)\\ R(l, u) &=& R(l, m)R(m, u) \end{eqnarray} \]

このとき、$m\simeq \left\lfloor\dfrac{l+u}{2}\right\rfloor$ とすることで上記のような計算量に落ち着く。

具体例の引用として $\arctan(1/x)$ の計算をマクローリン展開計算を取る。 \[ \arctan\frac{1}{x} \simeq \frac{P(0,N)}{Q(0,N)} =\sum_{k=0}^{N-1}\frac{(-1)^k}{(2k+1)x^{2k+1}} \] \[ \Rightarrow \left\{ \begin{eqnarray} P(l, u) &=& R(l, u) \sum_{k=l}^{u-1} (-1)^{k+1} \frac{x^{u-1-k}}{2k+3}\\ Q(l,u) &=& x^{2(u-l)} \prod_{k=l}^{u-1}(2k+3)\\ R(l,u) &=& \prod_{k=l}^{u-1}(2k+3) \end{eqnarray} \right. \]

詳細版

BS 法は上記の形がよく知られているが、本来はより広範囲に適用できるアルゴリズムである。 [FT05]

求めようとする級数が \[ S = \sum_{n=0}^{\infty} \frac{a(n)}{b(n)} \frac{p(0)\cdots p(n)}{q(0)\cdots q(n)} \quad {\rm または}\quad U = \sum_{n=0}^{\infty} \frac{a(n)}{b(n)} \left(\frac{c(0)}{d(0)}+\cdots+\frac{c(n)}{d(n)}\right) \frac{p(0)\cdots p(n)}{q(0)\cdots q(n)} \]

の形で表わされているとする。このとき、この部分級数 \[ S = \sum_{n_1 \lt n\lt n_2} \frac{a(n)}{b(n)} \frac{p(n_1)\cdots p(n)}{q(n_1)\cdots q(n)} \quad {\rm または}\quad U = \sum_{n_1 \lt n \lt n_2} \frac{a(n)}{b(n)} \left(\frac{c(n_1)}{d(n_1)}+\cdots+\frac{c(n)}{d(n)}\right) \frac{p(n_1)\cdots p(n)}{q(n_1)\cdots q(n)} \]

を直接計算するのではなく $P=p(n_1)\cdots p(n_2-1)$、$Q=q(n_1)\cdots q(n_2-1)$、 $B=b(n_1)\cdots b(n_2-1)$、$D=d(n_1)\cdots d(n_2-1)$、 $C=D(c(n_1)/d(n_1)+\cdots+c(n_2-1)/d(n_2-1))$ を考えて $T=BQS$ や $V=DBQU$ を求めてから $S=R/BQ$ や $U=V/DBQ$ を導く。 $m=\lfloor(n_1+n_2)/2\rfloor$ として $n_1\leqq n\lt m$ に属する項を $B_l$、$P_l$、$Q_l$、$T_l$、$D_l$、$C_l$、$V_l$、 逆に $m\leqq n\lt n_2$ に属する項を $B_r$、$P_r$、$Q_r$、$T_r$、$D_r$、$C_r$、$V_r$、 とするとき

\[ B = B_lB_r,\ P = P_lP_r,\ Q = Q_lQ_r,\ T = B_rQ_rT_l + B_lP_lT_r \] \[ D = D_lD_r,\ C = C_lD_r + C_rD_l,\ V = D_rB_rQ_rV_l + D_rC_lB_lP_lT_r + D_lB_lP_lV_r \]

で再帰的に計算する。 また以下のような C++ によるサンプルコード(コメント部分は日本語訳)も掲載されている。[FT05]

// S を計算する方
// ユーザーが初期値を入れるものとする
struct {
  bigint *a, *b, *p, *q;
} abpq_series;

// 分数の形で結果を保持する
struct {
  bigint P, Q, B, T;
} abpq_series_result;

// Binary Splitting による計算
void sum_abpq(abpq_series_result& r, int n1, int n2, const abpq_series& arg) {
  switch(n2 - n1) {  // 級数の長さで分類
  case 0:
    error_handler("summation devide", "n2-n1 should be > 0.");
    break;

  case 1: // the result at the point n1
    r.P = arg.p[n1];
    r.Q = arg.q[n1];
    r.B = arg.b[n1];
    r.T = arg.a[n1] * arg.p[n1];
    break;

  // case 2, 3, 4 は簡略化のため省略

  default: // 一般の場合
    abpq_series_result L, R;
    int nm = (n1 + n2) / 2;
    sum_abpq(L, n1, nm, arg);
    sum_abpq(R, nm, n2, arg);

    r.P = L.P * R.P;
    r.Q = L.Q * R.Q;
    r.B = L.B * R.B;
    r.T = R.B * R.Q * L.T + L.B * L.P * R.T;
    break;
  }
}
// U を計算する方
// ユーザーが初期値を入れるものとする
struct {
  bigint *a, *b, *p, *q, *c, *d;
} abpqcd_series;

// 分数の形で結果を保持する
struct {
  bigint P, Q, B, T, C, D, V;
} abpqcd_series_result;

// binary splitting による計算
void sum_abpqcd(abpqcd_series_result& r, int n1, int n2, const abpqcd_series& arg) {
  switch(n2 - n1) {  // 級数の長さで分類
  case 0:
    error_handler("summation devide", "n2-n1 should be > 0.");
    break;

  case 1: // the result at the point n1
    r.P = arg.p[n1];
    r.Q = arg.q[n1];
    r.B = arg.b[n1];
    r.T = arg.a[n1] * arg.p[n1];
    r.C = arg.c[n1];
    r.D = arg.d[n1];
    r.V = arg.a[n1] * arg.c[n1] * arg.p[n1];
    break;

  // case 2, 3, 4 は簡略化のため省略

  default: // 一般の場合
    abpqcd_series_result L, R;
    int nm = (n1 + n2) / 2;
    sum_abpqcd(L, n1, nm, arg);
    sum_abpqcd(R, nm, n2, arg);

    r.P = L.P * R.P;
    r.Q = L.Q * R.Q;
    r.B = L.B * R.B;
    bigint tmp = L.B * L.P * R.T;
    r.T = R.B * R.Q * L.T + tmp;
    r.D = L.D * R.D;
    r.C = L.C * R.D + R.C * L.D;
    r.V = R.B * (R.B * R.Q * L.V + L.C * tmp) + L.D * L.B * L.P * R.V;
    break;
  }
}

各関数に対応する項

上記の形で書かれても今一つ項が見えない…というより分かりづらいので [FT05] 中にある項を表にまとめてみた。 表の中に現れる変数 $x$ は有理数で $x=u/v$ と表されるものとする。 また、$p(n)$ や $q(n)$ は $n=0$ で特殊な値を取ることが多いので別項を立てたが、 一般項で記述可能な場合は - を入れている。

$f(x)$ $a(n)$ $b(n)$ $p(0)$ $q(0)$ $p(n)$ $q(n)$ 備考
$\exp(x)$ 1 1 1 1 $u$ $nv$
$\ln(1+x)$ 1 $n+1$ $u$ - $-u$ $v$
$\sin(x)$ 1 1 $u$ $v$ $-u^2$ $2n(2n+1)v^2$
$\cos(x)$ 1 1 1 1 $-u^2$ $2n(2n-1)v^2$
$\arctan(x)$ 1 $2n+1$ $x$ - $-x^2$ 1 何故か $u$、$v$ を使った形で与えられてない
$\dfrac{\sqrt{C}}{12\pi}$ $A+nB$ 1 1 1 $-(6n-5)(2n-1)(6n-1)$ $n^3C/24$ Chudnovskyの公式
各定数はリンク先と一緒。
カタラン数
$G$※1
1 $2n+1$ 1 1 $n$ $4n+2$ ※1
$2\zeta(3)$ ※2 $205n^2+250n+77$ 1 1 - $-n^5$ $32(2n+1)^5$
  • ※1 $G=\dfrac{3}{8}{\displaystyle \sum_{n=0}^{\infty}}\dfrac{1}{\left(2n \atop n\right)(2n+1)^2}+\dfrac{\pi}{8}\log(2+\sqrt{3})$ の $\sum$ 部分だけを求めている。