円周率.jp > 多数桁計算 > Binary splitting

Binary splitting

Binary splitting method

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

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

簡易説明版

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

計算しようとする対象が
\frac{P(0,N)}{Q(0,N)}\sum_{n=0}%5e{N-1} \frac{a(n)}{b(n)} \frac{p(0)\cdots p(n)}{q(0)\cdots q(n)}
のような形になるよう, 適切な式 P(lu), Q(lu), R(lu) を設定する. ([FT05]中でいくつか具体例が挙げられている)

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

P(lu) = P(lm) Q(mu) + R(lm) P(mu)
Q(lu) = Q(lm) Q(mu) ,   R(lu) = R(lm) R(mu)

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

arctan\frac{1}{x}\frac{P(0,N)}{Q(0,N)}\sum_{k=0}%5e{N-1} \frac{(-1)%5ek}{(2k+1)x%5e{2k+1}}
P(lu) = R(lu) \sum_{k=l}%5e{u-1} (-1)k+1 \frac{x%5e{u-1-k}}{2k+3}
Q(lu) = x2(u-l) \prod_{k=l}%5e{u-1}(2k+3) ,   R(lu) =\prod_{k=l}%5e{u-1}(2k+3)

初出論文通り

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

S\sum_{n=0}%5e{\infty} \frac{a(n)}{b(n)} \frac{p(0)\cdots p(n)}{q(0)\cdots q(n)}  または  U\sum_{n=0}%5e{\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<n<n_2} \frac{a(n)}{b(n)} \frac{p(n_1)\cdots p(n)}{q(n_1)\cdots q(n)}  または  U\sum_{n_1<n<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)}

を直接計算するのではなく Pp(n1)…p(n2−1), Qq(n1)…q(n2−1), Bb(n1)…b(n2−1), Dd(n1)…d(n2−1), CD(c(n1)/d(n1)+…+c(n2−1)/d(n2−1)) を考えて TBQSVDBQU を求めてから ST/BQUV/DBQ を導くようにする. m=[(n1n2)/2] として n1nm に属する項を BlPlQlTlDlClVl, 逆に mnn2 に属する項を BrPrQrTrDrCrVr とするとき

BBlBrPPlPrQQlQrTBrQrTlBlPlTr
DDlDrCClDrClDrVDrBrQrVlDrClBlPlTrDlBlPlVr

で段階的に計算できると書いている.また,以下の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は有理数で xu/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(x) 1 n+1 u - u v 自然対数,x−1=u/v
sin(x) 1 1 u v u2 2n(2n+1)v2
cos(x) 1 1 1 1 u2 2n(2n−1)v2
arctan(x) 1 2n+1 x - x2 1 何故かuvを使った形で
与えられてない
\frac{\sqrt{C}}{12\pi} AnB 1 1 1 −(6n−5)(2n−1)(6n−1) n3C/24 Chudnovskyの公式
各定数はリンク先と一緒.
カタラン数
G ※1
1 2n+1 1 1 n 4n+2 ※1
2ζ(3) ※2 205n2+250n+77 1 1 - n5 32(2n+1)5
  • ※1 G\frac{3}{8}\sum_{n=0}%5e{\infty}\frac{1}{\left(2n \atop n\right)(2n+1)%5e2}\frac{\pi}{8}log(2+\sqrt{3}) のΣ部分だけを求めている.