格子点の数え上げの高速化

凸関数・凹関数で囲まれる範囲の格子点の個数の数え上げなどを、凸包を用いて高速化するテクニックです。

半径 $n$ の円に含まれる格子点の個数を $O(n^{2/3})$ ほど、$1$ から $n$ までの約数の総和を $O(n^{1/3} \log{(n)})$ ほどで求めることができます。

考え方

格子点の個数は、floor 関数の和を用いて求められる場合が多いです。

例えば、半径 $n$ の円の第一象限内の格子点の個数は $\sum_{x=1}^{n} \left\lfloor\sqrt{n^2 - x^2}\right\rfloor$、$1$ から $n$ までの約数の個数の総和(つまり、$y = n/x$ の第一象限内の格子点の個数)は $2 \cdot \left(\sum_{i=1}^{\lfloor{\sqrt{n}}\rfloor} \left\lfloor\frac{n}{i}\right\rfloor\right) - \lfloor{\sqrt{n}}\rfloor^2$ などを用いて計算することが多く、それぞれ $O(n)$, $O(\sqrt{n})$ 時間で求まります。

ただ、これらの関数 ($f(x) = \sqrt{n^2 - x^2}$, $f(x) = n/x$) は傾きが緩やかに変化する $x$ の範囲(例えば、$|f''(x)| < 1$ の範囲)が比較的広く、整数点集合 $P(n) := \{ (x, \lfloor f(x) \rfloor) \mid 1 \le x \le n\}$ に対する凸包の端点数は少ないです。実際、$f(x) = \sqrt{n^2 - x^2}$ は $\Theta(n^{2/3})$ 個、$f(x) = n/x$ は $\Theta(n^{1/3} \log(n))$ (?) 個ほどの頂点によって凸包を構成できます。

凸包が取れると、元の問題は(関数が凸/凹なので)各台形中の格子点を数える問題となり、それらは(傾きが既約分数の場合)$O(1)$ で求められるので、凸包の端点数が少なければ計算量が改善します。

よって、凸包を $O($凸包の端点数$)$ ほどの計算量で構築できるとよいのですが、これは Stern–Brocot tree を用いて実現できます。

求め方

この方法の初出は Problem 540 - Project Euler での、Animus (+ Lucy_Hedgehog) さんの書き込みです。

例として、凹関数 $f(x) := \sqrt{n^2 - x^2}$ に対する整数点集合 $P(n) = \{ (x, \lfloor f(x) \rfloor) \mid 1 \le x \le n\}$ について凸包を高速に構築する方法を書きます。

今、凸包の端点 $(x_0, y_0)$ を見つけていると仮定し、次の端点 $(x_1, y_1)$ ($x_1 > x_0$) を見つけたいとします(正確には、線分 $(x_0, y_0)$, $(x_1, y_1)$ 上に存在する次の整数座標)。この時、$\frac{a}{b} < \frac{y_0 - y_1}{x_1 - x_0} \le \frac{c}{d}$ を満たす非負整数 $a, b, c, d$ ($bc - ad = 1$) も所持しているとします。

このとき、次のようにすると Stern-Brocot tree の性質により $(x_1, y_1)$ を求められます。

  1. $p := a + c, q := b + d$ とし、
    • $(x_0 + q)^2 + (y_0 - p)^2 \le n^2$ の場合、$c \gets p$, $d \gets q$ として、ステップ 1. に戻ります。
    • $(x_0 + q)^2 + (y_0 - p)^2 > n^2$ の場合、
      • $-f'(x_0 + q) \ge \frac{c}{d}$ を満たしていると、今後 $(x_0 + q)^2 + (y_0 - p)^2 \le n^2$ となることはないので、$(x_1, y_1) := (x_0 + d, y_0 - c)$ として処理を打ち切ります。
      • それ以外の場合、$a \gets p$, $b \gets q$ として、ステップ 1. に戻ります。

実際には、$a, b, c, d$ を効率よく見つけるために Stern-Brocot tree 上の遷移履歴を stack に積んでおく必要もあり、例えば以下のように実装できます。

計算量は(少なくとも)$\Omega(n^{2/3})$ かかりますが、実用上 $\Theta(n^{2/3})$ と見なしても問題なさそうです。

(しっかりとした証明をできていないので、以下 $\Omega$ で書きます。)

例題

  • DIVCNT1: http://www.spoj.com/problems/DIVCNT1/
    • $\sigma_0(n)$ の総和に対応します。$f''(x) = 2n/x^3$ なので、$\sqrt[3]{2n} < x \le \sqrt{n}$ の部分について、凸包を取ります。
  • 他にも 60, 90, 120 度の角を持つ整数辺三角形の数え上げなども、格子点の数え上げ問題に帰着できます。

応用

前の高速化を使うと、$$S_k(n) := \sum_{i=1}^n \sigma_k(i) = \sum_{i=1}^n i^k \left\lfloor\frac{n}{i}\right\rfloor$$ が $\Omega(k^2 n^{1/3} \log{n})$ 回の演算で求まります。

格子点の個数を数える際(例えば、$S_{0} (n)$)では、stack には $(a, b)$ のみを積んでいくのですが、$S_{k} (n)$ を求める際には $(a, b, T_1(a, b), \dots, T_k(a, b))$ を stack に積みます。ここで、$T_k(a, b) := \sum_{x=0}^{b-1}x^k\left\lfloor\frac{ax}{b}\right\rfloor$ です。また、$P_k(b) := \sum_{x=0}^{b-1} x^{k}$ と定義します。

これは、$S_k(n)$ を凸包を用いて台形に分割して求める場合、各台形(傾き $a/b$)での部分和は $T_0(a, b), T_1(a, b), \dots, T_k(a, b)$ と $P_0(b), P_1(b), \dots, P_k(b)$ の値を用いて $O(k)$ 回の演算で計算できるためです(前計算ありの状態で、$P_i(b)$ をその場で計算すると合計 $O(k^2)$)。

$(a, b, T_1(a, b), \dots, T_k(a, b))$ と $(c, d, T_1(c, d),\dots, T_k(c, d))$ から $(a + c, b + d, T_1(a + c, b + d), \dots, T_k(a + c, b + d))$ を得る際には、前の記事の方法を用いることができ、(愚直な畳み込みにより)$O(k^2)$ 回の演算で済みます。

よって、$S_k(n)$ は $\Omega(k^2 n^{1/3} \log{n})$ 回の演算で求めることができます。

例題