テク名が長すぎませんか？
これはなに
https://min-25.hatenablog.com/entry/2018/05/03/145505min-25.hatenablog.com
を書いてみましたエッセイです(??) 上の記事を読んでいなくても話がわかるように書きます
モチベーション
例えば、次のような問題を考えてみましょう。
自然数 $x$ について、$x$ の正の約数の個数を $\sigma(x)$ と書くことにします。
自然数 $N,M(1\leq N\leq M)$が与えられるので、次の値を求めてください。
この問題は、$f(X)=\sum_{i=1}^X\sigma(i)$ に対して $f(M)-f(N-1)$ を求めることで計算でき、 $f(X)$ は（主客転倒と呼ばれがちなテクニックを用いて）次のように書けます。
これは、凸な領域 $S_X:=\{(x,y)\in\mathbb{R}^2|0\lt x\leq X,0\lt y\leq X,xy>X\}$ に含まれる格子点 $S_X\cap\mathbb{Z}^2$ の要素数を数えられれば解けます。
この格子点はちょうどよく凸っぽい性質を残しているので、もしこの格子点たちの凸包を生成する点だけがすべてわかっているなら、$O( (\text{凸包の頂点数}) )$ や $O((\text{凸包の頂点数})\log X)$ くらいの時間で点の個数を数え上げることができます（このパートは下の問題とほぼ同じです）。
$f(X)$ は普通に計算すると $\Theta(X)$ や $\Theta(\sqrt{X})$ 時間がちですが、凸包の頂点数はそれより小さくて $O(\sqrt[3]{X}\log X)$ 程度な1ので、真に高速です。
速いとうれしいので、やっていきましょう
凸集合の中の格子点の凸包の頂点の性質
まず、凸包について考えれば十分であることを示します。
$C\subseteq\mathbb{R}^2$ が凸集合であるとします。
$C$ 内の格子点全体 $C_\mathbb{Z}=C\cap\mathbb{Z}^2$ が「凸っぽい」性質をもっていること、つまり「$C_\mathbb{Z}$ は、その凸包に含まれる格子点全体である」ことを示します。
$C_\mathbb{Z}$ の要素が凸包に含まれないことはないので、凸包に含まれる格子点であって $C_\mathbb{Z}$ の要素でないものが存在しないことをいいます。
$C_\mathbb{Z}$ の凸包は $C_\mathbb{Z}$ を含む凸集合全体の共通部分です。 $C$ は $C_\mathbb{Z}$ を含む凸集合なので、$C$ は $C_\mathbb{Z}$ の凸包を完全に含み、凸包に含まれる格子点は $C$ に含まれます。 示されました。
なので、$C$ に含まれる格子点について知りたい場合は $C$ 内の格子点全体の凸包について知れば十分であることがわかりました。 逆に、$C_\mathbb{Z}$ の凸包について知りたいとき、$C$ についての情報があれば大丈夫そうということもわかりました。
本編
凸集合 $C$ が与えられたとき、$C_\mathbb{Z}=C\cap\mathbb{Z}^2$ の凸包の頂点を列挙することを考えます。
ここでは、「いちばん上」から「いちばん右」へ向かう部分のみについて考え、$(0,y)$ から $(x,0)$ まで移動していくことを考えます。 適宜平行移動したり反転したりしておいてください。
主に次の2つのパートからなります。
- 今の傾きのまま外に出ないギリギリまで進む
- 次の傾きを計算する
1は簡単（にぶたんで $O(\log X)$ くらい、閉じた式で書けるならそれを計算するのにかかる時間くらい）なので、2について詳しく書きます
今いる点を $(x,y)$ として、$(x+a,y-b)\in C,(x+c,y-d)\notin C,bc-ad=1$ なる $a,b,c,d$ がわかっているとします。 今回の状況では例えば $(a,b,c,d)=(0,1,1,0)$ が常にこれを（最後に点が $(x,0)$ に到達したとき以外）満たします。 これは Stern-Brocot Tree の部分木を探索することに対応しています。
求めるもの（次の傾き）は、$(x+\alpha,y-\beta)\in C$ を満たす互いに素な非負整数の組 $(\alpha,\beta)$ であって $\dfrac{\alpha}{\beta}$ が最大となるもの2です。
$(a,b)$ が求めるものである場合と、そうでない場合があります。 $(a,b)$ が求めるものであるとは、傾きが $(a,b)$ と $(c,d)$ の間にあるような格子点はすべて $C$ に含まれないということです。
これから先に見ることになる格子点 $(x+\alpha,y-\beta)$ は、すべて $\dfrac{a}{b}\lt\dfrac{\alpha}{\beta}\lt\dfrac{c}{d}$ を満たします。
このような格子点の存在する範囲を厳しめに評価することで、$C$ に含まれるものが存在しないことを確かめられます。
Stern-Brocot Tree の性質から、$\dfrac{a}{b}$ と $\dfrac{c}{d}$ の間にある既約分数はすべてある互いに素な正整数 $(n_1,n_2)$ に対して $\dfrac{n_1a+n_2c}{n_1b+n_2d}$ と書けます。 なので、これから先に見る格子点はすべて $(x+a+c+\alpha,y-b-d+\beta)\ (\alpha=n_1^\prime a+n_2^\prime c,\beta=n_1^\prime b+n_2^\prime d,0\leq n_1^\prime,n_2^\prime)$ の形で書けることがわかります。
よって、$(x+(t+1)a+c,y-(t+1)b-d)\ (0\leq t)$ が $C$ と共通点を持たないとき、$(a,b)$ が求める次の傾きです。
そうでないときは、Stern-Brocot Tree の子に動けばいいです。 $(x+a+c,y-b-d)$ が $C$ に含まれるかどうかで場合分けを行います。
この繰り返しは、求める傾きの Stern-Brocot Tree における深さの回数まで行われます（例えば、次向かうべき向きが $(13,5)$ なら下の図のように辿り、$(a,b,c,d)=(0,1,1,0)\to(1,1,1,0)\to(2,1,1,0)\to(2,1,3,1)\to(5,2,3,1)\to(5,2,8,3)\to(13,5,8,3)$ もしくは $(5,2,13,5)$ のように7回の繰り返しを行います）。
出てくる傾きの Stern-Brocot Tree における深さの最大値は $x^2+y^2\leq n$ では $\Omega(\sqrt[4]{n})$ 、$xy\geq n$ では $\Omega(n)$ となります。
これらの場合でも Stern-Brocot Tree の同じ向きの遷移をまとめることで遷移の回数が $O(\log n)$ で抑えられるようになり、いい感じです。
これを繰り返して探索をしていきます。 一度探索をした後の次の探索で求める傾きはさっき求めた傾きと近いことが予想されるので、Stern-Brocot Tree に潜ってきた道を残しておくことで高速化になりそうです（ここの計算量の証明をしていません 毎回一から探索しても $O(\log n)$ 倍になるだけなのできっと大丈夫です）。
実装方針/作例
ざっくり次の3つの機能があれば動くことがわかります。
- 点 $(x,y)$ が領域の中か外かを判定する
- 点 $(x,y)$ から $(a,-b)$ 方向（$(x+ta,y-tb)$ です）に伸びる半直線と境界とが交わるか判定する
- 点 $(x,y)$ から $(a,-b)$ 方向に伸びる半直線と境界とが交わることがわかっているとき、境界上にある $(x+ta,y-tb)$ のうちもっとも小さい $t\geq0$ を求める
これらの関数 is_inside、is_cross、next_in/next_out （始点 $(x,y)$ が領域の中にあるかで関数を分けたほうが気持ちがよかったのでそうしています $t$ の丸め方の方向も違うので） を使って、一回分のロジックはこんな風に書くことができます。
now_point には今見ている点（凸包の頂点）、stern_brocot_tree には探索の履歴、result には今まで見つかった凸包の頂点が入っています。
auto& [x, y]{now_point}; // subtree between a/b and c/d auto [a, b, c, d]{stern_brocot_tree.back()}; stern_brocot_tree.pop_back(); // next vertex of convex hull const auto n_upper{next_out(x, y, a, b, N)}; x += n_upper * a; y -= n_upper * b; result.emplace_back(now_point); // backtrack while(!empty(stern_brocot_tree)){ tie(a, b, c, d) = stern_brocot_tree.back(); if(is_inside(x + a, y - b, N))break; stern_brocot_tree.pop_back(); } if(empty(stern_brocot_tree))break; // search forward while(y >= d){ if(y >= b + d && is_inside(x + a + c, y - b - d, N)){ const auto n_lower{next_out(x + a + c, y - b - d, c, d, N)}; tie(a, b, c, d) = make_tuple(a + c * (n_lower + 1), b + d * (n_lower + 1), a + c * (n_lower + 2), b + d * (n_lower + 2)); }else if(is_cross(x + c, y - d, a, b, N)){ const auto n_lower{next_in(x + c, y - d, a, b, N)}; if(y < b * n_lower + d || !is_inside(x + a * n_lower + c, y - b * n_lower - d, N))break; tie(a, b, c, d) = make_tuple(a * n_lower + c, b * n_lower + d, a * (n_lower - 1) + c, b * (n_lower - 1) + d); }else break; stern_brocot_tree.emplace_back(a, b, c, d); }
subtree を特定するのには $\dfrac{a+c}{b+d}$ があれば十分だと思うんですけど、$\dfrac{a+c}{b+d}$ から $\dfrac{a}{b},\dfrac{c}{d}$ を復元する高速な方法がわからなかったのでそのままにしています 情報求ム
あとは、これに合うように初期値と各関数を定めていけばいいです。
例題
求める自然数の組 $(x,y)$ は必ず $x^2+y^2\leq N$ の中の格子点全体の凸包の頂点になっています（背理法で示せます）。 なので、凸包の頂点のみを見れば十分で、今回のテクが使えます。
考えるべき領域が四分円 $x^2+y^2\leq N,x\geq0,y\geq0$ なので、is_inside、is_cross、next_in/next_out は次のような判定をすればいいです。
- $x^2+y^2\leq N\wedge x\geq0\wedge y\geq0$
- $\dfrac{|ax-by|^2}{a^2+b^2}\leq N$
- $\dfrac{by-ax\pm\sqrt{(a^2+b^2)N-(ay+bx)^2}}{a^2+b^2}$
大きな $n$ について計算するとき、オーバーフローに気をつけてください。
点の初期値は $(x,y)=(0,\left\lfloor\sqrt{N}\right\rfloor)$ 、stern_brocot_tree には Stern-Brocot Tree 全体にあたる $(0,1,1,0)$ などを入れておけばいいです。
初期値を $(x,y)=\left(\sqrt{\dfrac{N}{2}},\sqrt{\dfrac{N}{2}}\right)$ 付近にして stern_brocot_tree を $(0,1,1,1)$ から始めると2倍の高速化になります。
これを実装すると、次のようになります。
だいたい $2^{64}$ くらいの数 $18446744072107595890$ を入れると 380 ms くらいかかって結果が返ってきます。
上で述べた定数倍高速化をするとしっかり2倍で 190 ms くらいになります。
シンプルな $\Theta(\sqrt{N})$ を書くとコードテストではタイムアウト、手元では2分ほどかけて答えを出力しました。
はやくなってうれしいね 愚直を書くのが上手な人にコードテスト上で通されると困った気持ちにはなります
例題 2
自然数 $x$ について、$x$ の正の約数の個数を $\sigma(x)$ と書くことにします。
自然数 $N,M(1\leq N\leq M)$が与えられるので、次の値を求めてください。
一番最初に出した問題です
$f(X)=\sum_{i=1}^X\sigma(i)$ を求めたいとき、考える領域は $xy\gt X,x\gt0,y\gt0$ です。 今は格子点だけ注目しているので $xy\geq X+1,x\gt0,y\gt0$ と書き換えていいです。 $X+1$ を改めて $N$ とすることにしましょう。
これは下に凸な関数なのでちょっとだけ変えないとうまく回りませんが、ちょっと変えるだけでうまく回ります。 具体的には、部分木を表す $(a,b,c,d)$ を $(c,d,a,b)$ のように保管するようにするとほぼ同じコードでうまく動くようになります。
初期点を $(1,N)$、stern_brocot_tree を $(1,0,0,1)$ から始めることでうまくやっていくことができます。
is_inside、is_cross、next_in/next_out は次のような判定をすればいいです。
- $xy\geq N\wedge x\geq0\wedge y\geq0$
- $ab=0$ のとき $bx-ay\gt0$、$ab\neq0$ のとき$(bx+ay)^2\leq 4abN$
- $a=0$ のとき $\dfrac{xy-N}{bx}$、$b=0$ のとき $\dfrac{N-xy}{ay}$、$ab\neq0$ のとき $\dfrac{ay-bx\pm\sqrt{(ay+bx)^2-4abN}}{2ab}$
こちらも大きい $N$ で計算する場合はオーバーフローに気をつけましょう。
こっちにも対称性があるので、$\left(\sqrt{N},\sqrt{N}\right)$ あたりから始めることで定数倍が高速になります。
実装は次のようになります。
判定パートをオーバーフローさせない方法が思いつかなかったので unsigned long と __uint128_t 固定になっています。
また、凸包を構築する途中に「注目している辺の上にいくつ格子点があるか」がわかっているので、構築しながら答えも求めると $\gcd$ 一回分の計算量が浮いてお得です。
答えを unsigned long で計算しているので適当な $N,M$ を放り込むとオーバーフローすると思います。←??
1つめの例題ほど速くはなくて、$2^{64}$ くらいの数を入れるとコードテストで間に合わなかったりします（数が大きいところでバグっているのかもしれません(??)）。 だいたい $10^{18}$ 前後の数 $N=662891320790343374,M=1808274646968721915$ を入れると 3秒ほどで答えが出ます。
__uint128_t などを経由していて定数倍も見るからに重そうで、凸包の頂点数の評価もひとつめの $O(\sqrt\[3\${n})] から $O(\sqrt\[3\${n}\log n)] になっていて、実際重いです。
これも定数倍高速化はきっちり 2 倍速になってだいたい 1.5 秒で答えが出ます。
例題 3
例題 2 の領域で例題 1 みたいなことをすればいいです。
いかがでしたか？
🦑←これはイカ
オーバーフロー回避策、高速化テクなどあればどしどしお寄せください 宛先はありません
-
https://www.ams.org/journals/tran/1963-106-02/S0002-9947-1963-0143105-7/S0002-9947-1963-0143105-7.pdf の主定理を $n=2$ 次元の場合について使うと面積が $V$ の狭義凸集合の境界に格子点が $N$ 点存在するとき $N=O(\sqrt[3]{V})$ が言えるので、全体をうまく$O(\log X)$ 個の部分に分割して $O(\sqrt[3]{X}\log X)$ が示せます（thanks to Min_25 さん）
↩x の範囲を (sqrt(n)/2^{k+1}, sqrt(n)/2^k] (k = 0, 1, ...) で分けた後、https://t.co/2Cykibtnha で O(n^{1/3} log(n)) までは示せるかな (?) と思っています（Θ かどうかはわからないですが）
— Min_25 (@min_25_) 2021年7月13日 - ただし、$\beta=0$ のとき、そうでないどのような組よりも大きいとします↩
