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

Binary splitting method

日本語で言うなら二分割法，とでもいうのだろうか． [FT05]によると各関数や定数についていえば Karatsuba が同様の概念を提案していたらしいが， [FT05]に一般化した2種類の形で纏められている．

分割統治を行うことで，n 桁乗算の計算量を M(n) とするときの計算量が O(logn M(n logn)) となる．[FT06]

[FT06] では手法が簡略に示されているのでそちらの形で紹介する．

計算しようとする対象が

のような形になるよう， 適切な式 P(l，u)， Q(l，u)， R(l，u) を設定する． （[FT05]中でいくつか具体例が挙げられている）

それらの具体的計算については， $m=\left\lfloor\frac{l+u}{2}\right\rfloor$ として以下のような再帰的な計算ができるようにする．

具体例の引用として $\arctan(1/x)$ の計算をマクローリン展開計算を取る． $l$ と $u$ が等しい場合は定義式に基づいた計算を行えば具体的な計算ができる．

まずは[FT05]に見られる binary splitting の方法を見てみる． 求めようとする級数が

の形で表わされているとする．このとき，この部分級数

を直接計算するのではなく P ＝ p(n₁)…p(n₂−1)， Q ＝ q(n₁)…q(n₂−1)， B ＝ b(n₁)…b(n₂−1)， D ＝ d(n₁)…d(n₂−1)， C ＝ D(c(n₁)/d(n₁)＋…＋c(n₂−1)/d(n₂−1)) を考えて T＝BQS や V＝DBQU を求めてから S＝T/BQ や U＝V/DBQ を導くようにする． m＝[(n₁＋n₂)/2] として n₁≦n＜m に属する項を B_l，P_l，Q_l，T_l，D_l，C_l，V_l， 逆に m≦n＜n₂ に属する項を B_r，P_r，Q_r，T_r，D_r，C_r，V_r とするとき

で段階的に計算できると書いている．また，以下のC++によるサンプルコード（コメント部分は日本語訳しました）も掲載されている．

// 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で特殊な値を取ることが多いので別項を立てたが，一般項で記述可能な場合は-を入れている．

• ※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$ 部分だけを求めている．