2022-03-06

包除原理の基礎的な話

競技プログラミング

はじめに

包除原理を知っていますか？私は知っています

知ってはいますが, 全然使いこなせていないし使うたびに頭が壊れてしまって「答えが合ったので, ﾖｼ!!」と提出してしまうので, ちゃんと勉強しようと思いました

この記事はそのまとめです

本文

とりあえずググりましょう

一番上に wikipedia 先生が出てきたので開いてみると, このような式が乗っています

$\left\vert \displaystyle \bigcup _ {i =1} ^ n A _ i \right\vert = \displaystyle\sum _ {k=1} ^ n (-1) ^ {k - 1} \displaystyle\sum _ { \vert J \vert = k} \left\vert \displaystyle\bigcap _ {j \in J} A _ j \right\vert$

複雑そうに見えるかもしれませんが, 「積集合をいい感じに足し引きすると和集合になる」というだけです

証明はいつもお世話になっております高校数学の美しい物語さんを見るといいでしょう

この式を見ればわかる通り, これは積集合から和集合を求めています

つまり, ある集合 $S$ を決め打ったときに「少なくとも $S$ を含む場合」が簡単にわかるときに「少なくとも一つを含む場合」を求めようとしています

また, $A _ 1, A _ 2, ... , A _ n$ の上位集合を $X$ とすると

$X - \left\vert \displaystyle \bigcup _ {i =1} ^ n A _ i \right\vert = X - \displaystyle\sum _ {k=1} ^ n (-1) ^ {k - 1} \displaystyle\sum _ { \vert J \vert = k} \left\vert \displaystyle\bigcap _ {j \in J} A _ j \right\vert \ = \displaystyle\sum _ {k=0} ^ n (-1) ^ {k} \displaystyle\sum _ { \vert J \vert = k} \left\vert \displaystyle\bigcap _ {j \in J} A _ j \right\vert$

と書けます ( $0$ 個の積集合 = 上位集合全体と考えています)

これは元の式の余事象なので, 「一つも含まない場合」を意味します

ここまで集合の形で書きましたが, 競プロではだいたいこのままの形では出ません

競プロの問題文を要約すると, たいていの場合

いくつかの条件が与えられるこれらをすべて満たすようなものが何通りあるか求めよ

みたいな感じになります

このせいで, 「なんか包除するのかな～」という気持ちになっても「条件のうち, これとこれは満たすものを足し合わせて... ？」となって混乱してしまいます (僕だけですか？)

上に書いた通り, 包除原理で求めることができるのは和集合です

「すべて満たすようなもの」というのは積集合なのでこのままでは求められません

そのため, 逆に「条件を満たさないもの」を考えます

つまり,

いくつかの条件が与えられるこれらに一つも違反していないものは何通りあるか求めよ

という感じに言い換えるわけです

すると, この「これらに一つも違反していないもの」は「少なくともこれとこれは違反するような方法は何通りあるか」のようなものが求められれば包除原理で求められるわけです

おわりに

個人的に, この「条件を満たさないものを考える」というのがとても非直感的だなと感じてしまうため, いつも思考が妨げられていました

こうやってちゃんと考えて少しずつ慣れていけたらなと思います

参考文献

包除原理 - Wikipedia

包除原理の２通りの証明 | 高校数学の美しい物語

2021-10-01

Cipolla's algorithm (素数 mod 上の平方根)

競技プログラミング

はじめに

素数 mod 上の平方根を求める方法を勉強したので, 軽くまとめます

体とか同型とかそういう言葉を使わないように書いていきたいと思っています

ちゃんと勉強したい人は 37zigenさんのブログを読んでください

やりたいこと

素数 $p$ と整数 $a$ が与えられるので

$x ^ 2 \equiv a \mod p$

となるような $x$ を一つ求めてください

存在しないならそのことを報告してください

やりかた

めんどくさいのでここから下では $\mod p$ をつけません (ごめんなさい)

何も言わずに $\equiv$ って書いてあったら $\mod p$ だと思ってください

とりあえず

まずは自明なものをはじきましょう

というのも, 今回使う道具は奇素数にしか対応してないです

$p = 2$ のとき $x, a$ は $0$ か $1$ のどっちかですが,

$0 ^ 2 = 0, 1 ^ 2 = 1$ なので, $x = a$ です簡単ですね

ついでに $a = 0$ のときも先にはじきます $x = 0$ ですね

ということでここから先は $p$ は奇素数, $1 \le a \lt p$ とします

平方剰余

そしたらまずは, そもそも平方根が存在するのかを判定しましょう

まず用語の説明ですが,

「 $x ^ 2 \equiv a$ となる整数 $x$ が存在する」とき, 「 $a$ は法 $p$ で 平方剰余 である」といいます

オイラーの規準

さて, この判定をすることができる道具として, オイラーの規準 というものがあります

これは,

$a ^ {(p-1)/2} \equiv 1$ のとき $a$ は法 $p$ で平方剰余である

$a ^ {(p-1)/2} \equiv -1$ のとき $a$ は法 $p$ で平方剰余ではない

というものですこれは逆も成り立ちます

$\dfrac{p - 1}{2}$ がついてますねこれのために $p = 2$ をはじいたわけです

さて, $a ^ {(p-1)/2}$ が $1$ でも $-1$ でもない値になることはないのか？と思うかもしれません (僕は思いました)

これは, ならないことがわかります

フェルマーの小定理より, $p$ が素数で $a$ が $p$ と互いに素であるとき

$a ^ {p-1} \equiv 1$ となります

ということは,

$\left( a ^ {(p-1)/2} \right) ^ 2 \equiv a ^ {p-1} \equiv 1$

ですね

二乗して $1$ になる整数は $1$ と $-1$ しかないので, $a ^ {(p-1)/2}$ はそのどちらかになるということです

証明

$1$ か $-1$ になることはわかりましたが, なぜそれで判定できるのでしょう

まず, 「 $a$ は法 $p$ で平方剰余である」 $\Rightarrow a ^ {(p-1)/2} \equiv 1$ を示します

$a$ は法 $p$ で平方剰余なので, とある $x$ が存在して $x ^ 2 \equiv a$ です

ここで, $x$ も $p$ と互いに素なはずです (じゃないと $0$ になってしまうので)

ということで, フェルマーの小定理より $x ^ {p-1} \equiv 1$ です

$p-1$ は偶数なので $x ^ {p-1} \equiv (x ^ 2) ^ {(p-1)/2} \equiv a ^ {(p-1)/2}$ と書けますね

ということで $a ^ {(p-1)/2} \equiv 1$ です

次に, $a ^ {(p-1)/2} \equiv 1 \Rightarrow$ 「 $a$ は法 $p$ で平方剰余である」を示します

ここで, $r$ を $p$ に対する原始根とします

$r$ は原始根なのである自然数 $k$ が存在して $a \equiv r ^ k$ と書けるはずです

よって $a ^ {(p-1)/2} \equiv r ^ {k(p-1)/2}$ となります

ここで $a ^ {(p-1)/2} \equiv 1$ とすると, $r ^ {k(p-1)/2} \equiv 1$ です

$r$ は原始根なので, $k(p-1)/2$ は $p-1$ の倍数である必要があります

要らない気がしますがこれも証明してみると

$k(p-1)/2$ が $p-1$ の倍数でない場合,

$k(p-1)/2 = (p-1)A + B \,$ ( $A , B$ は非負整数, $0 \lt B \lt (p-1)$ ) と書けます

ということは

$r ^ {k(p-1)/2} = r ^ {(p-1)A + B} = r ^ {(p-1)A} \cdot r ^ B$

となるので

$r ^ {(p-1)A} \cdot r ^ B \equiv 1$ です

ここで, $r ^ {p-1} \equiv 1$ なので $r ^ B \equiv 1$ であるはずですがこれは $r$ が原始根であることに矛盾します

ということで $k(p-1)/2$ は $p-1$ の倍数です

つまり $k$ が偶数であるわけですが, $k = 2n$ としてみると

$a \equiv r ^ k \equiv r ^ {2n} \equiv (r ^ n) ^ 2$

と書けるので $a$ の平方根 $r ^ n$ が存在します

証明終わり！

Cipolla's algorithm

ということで解が存在するかどうかの判定はできました

ここからは存在するものとします

ここで, いい感じの整数 $b$ を使って

$f(x) = x ^ 2 - 2bx + a$

という方程式を作ってみます

$f(x) = 0$ の解は, 普通に解の公式を使うと

$x = b \pm \sqrt{b ^ 2 - a}$

と書けます

さて, $b ^ 2 - a$ が平方剰余である場合,

ある非負整数 $t$ が存在して $t ^ 2 \equiv b ^ 2 - a$ と書けるので $x \equiv b \pm t$ です

だからなんですか？

何もありません

何もないと困るので, $b ^ 2 - a$ は平方剰余じゃないとします

それはそれで困りました平方剰余じゃないってことは $\sqrt{b ^ 2 - a}$ は整数じゃないです

仕方ないのでそのまま残して計算することにします

どういうことかというと, $a + b\sqrt{X}, c + d\sqrt{X}$ とあったとすると

$(a + b\sqrt{X}) + (c + d\sqrt{X}) = (a + c) + (b + d)\sqrt{X}$

$(a + b\sqrt{X})(c + d\sqrt{X}) = (ac + bdX) + (ad + bc)\sqrt{X}$

という感じです

また, 項ごとに mod を取れることにします

$7 + 6\sqrt{X} \equiv 3 + 2\sqrt{X} \mod 4$

ということで計算していきます

書くのがめんどくさいので

$b + \sqrt{b ^ 2 - a} , b - \sqrt{b ^ 2 - a} = \alpha, \beta$ とおきます

さて, ここでこいつらの $p$ 乗を考えてみます

さっきまでならフェルマーの小定理から $a ^ p \equiv a$ みたいなことが言えましたが, もうそんなものは成り立ちません

$\alpha ^ p$ を二項定理で頑張ってみましょう

めんどくさいので $\sqrt{b ^ 2 - a} = X$ とおいて

$\alpha ^ p = (b + X) ^ p = \displaystyle\sum _ {k = 0} ^ p \binom{p}{k} b ^ k X ^ {p-k}$ と書けます

$\binom{p}{k}$ について考えてみます

定義より, $\binom{p}{k} = \dfrac{p!}{k! (p-k)!}$ と書けます

$1 \le k \le p-1$ としてみると, 分母の $k! (p-k)!$ に $p$ が登場しません

よって分子の $p!$ の $p$ が約分されないため, $\binom{p}{k}$ は $p$ の倍数になります

ということで mod を取ると $k = 0, p$ だけが残って

$\alpha ^ p \equiv b ^ p + X ^ p$ となります (あれ, 一年生の夢が現実になりました)

ここで, $b$ は普通に整数だったのでフェルマーの小定理から $b ^ p \equiv b$ です

$X$ の方は, $p$ が奇数なので $X ^ p = (X ^ 2) ^ {(p-1)/2} \cdot X$ と書けます

$X$ を $\sqrt{b ^ 2 - a}$ に戻すと

$(\sqrt{b ^ 2 - a}) ^ p = (b ^ 2 - a) ^ {(p-1)/2} \sqrt{b ^ 2 - a}$ となります

$b ^ 2 - a$ は平方剰余ではありませんでしたよね？

ここで, オイラーの規準を思い出してみましょう

$b ^ 2 - a$ は平方剰余ではないので $(b ^ 2 - a) ^ {(p-1)/2} \equiv -1$ と書けます

これを代入すると

$(\sqrt{b ^ 2 - a}) ^ p \equiv -\sqrt{b ^ 2 - a}$

と書けましたすっきりしましたね

そんなわけで, $\alpha ^ p \equiv b ^ p + X ^ p \equiv b - \sqrt{b ^ 2 - a} \equiv \beta$

となりました

ここで, 解と係数の関係を考えてみると, $\alpha \beta = a$ です

$\alpha ^ p \equiv \beta$ を代入すると $\alpha ^ {p + 1} \equiv a$ となります

ということで,

$(\alpha ^ {(p + 1)/2}) ^ 2 \equiv a$ です

平方根が見つかりました！！！！

いい感じに終わった風ですが, 問題があります

ちょっと前のところを思い出してみましょう

ここで, いい感じの整数 $b$ を使って

$f(x) = x ^ 2 - 2bx + a$

という方程式を作ってみます

$b$ ってどうやって見つければいいんだ...

さて, この $b$ に求められる条件というのは

$b ^ 2 - a$ が平方剰余にならない

ということだけです

実は, $a$ が平方剰余であるとして, $b$ をランダムに取ってきたときに $b ^ 2 - a$ が平方剰余にならない確率は $\dfrac{p-1}{2p}$ となります (これの証明はよくわかってません...)

追記: わかりました

$b ^ 2 - a$ が平方剰余になるような $b$ がいくつあるかを数えれば, その個数を $k$ 個とすると求めたい確率は $1 - \dfrac{k}{p}$ となります ( $b$ の候補は $p$ 個あるので)

さて, 今 $a$ は平方剰余であり $0$ でないので, とある整数 $m ( \not\equiv 0)$ を用いて $a \equiv m ^ 2$ と書けます

つまり $b ^ 2 - m ^ 2$ が平方剰余であるような $b$ の個数を数えればいいのですが, $m \not\equiv 0$ なので逆元が存在して, $m ^ 2 ( (m ^ {-1} b) ^ 2 - 1 )$ と書けます

$m ^ {-1} b = c$ とおくと, $c ^ 2 - 1$ が平方剰余かどうかを判定すればいいことになります

(オイラーの規準を思い出してあげるとわかりやすいのですが, 平方剰余である数 $\times$ 平方剰余である数は平方剰余である数, 平方剰余である数 $\times$ 平方剰余でない数は平方剰余でない数となります)

ここで唐突に $g(x) \equiv x ^ 2 - 2cx + 1$ という関数を考えてみます $g(x) \equiv 0$ の解が存在するためには $c ^ 2 - 1$ が平方剰余である必要があります (上の方でやったのとだいたい同じです) これは逆も成立するので同値です

整数 $t$ が解の一つ, つまり $g(t) \equiv 0$ であるとすると, $1$ が残ってしまうので $t \not\equiv 0$ であるといえますよってこいつにも逆元が存在して $c \equiv 2 ^ {-1} (t + t ^ {-1})$ と書けます

ということで $t + t ^ {-1}$ の値が何通りあるのかと変形できました

ここで $t + t ^ {-1} \equiv s + s ^ {-1}$ とすると, $t + t ^ {-1} - s - s ^ {-1} \equiv 0$ なのでこれを因数分解して $(t - s)(1 - (ts) ^ {-1}) \equiv 0$ , つまり $(t - s)(1 - ts) \equiv 0$ と書けます

つまり $t \equiv s$ または $t ^ {-1} \equiv s$ です

したがって, $t + t ^ {-1} \equiv s + s ^ {-1}$ なら $\{ t, t ^ {-1}\} = \{ s, s ^ {-1}\}$ となり, $\{ t, t ^ {-1}\}$ の個数が $k$ と一致します

$t$ は $p - 1$ 通り存在しますが, それらは $\pm 1$ を除いて $t \not\equiv s$ で $\{ t, t ^ {-1}\} = \{ s, s ^ {-1}\}$ となるペアが存在します

ということで $k = 2 + \dfrac{p - 3}{2} = \dfrac{p + 1}{2}$ となります

追記終わり

ということで $\dfrac{p-1}{2p}$ となりますが,

これはだいたい $\dfrac{1}{2}$ となります

ということで $n$ 回連続で失敗するような確率はだいたい $\dfrac{1}{2 ^ n}$ くらいと言えます

何が言いたいかというと, すぐ見つかるということです

まとめ

$p = 2$ や $a = 0$ の場合は early return
$a ^ {(p-1)/2}$ を計算して解が存在するか判定
$(b ^ 2 - a) ^ {(p-1)/2} \equiv -1$ になるような $b$ を頑張って探す
$(b + \sqrt{b ^ 2 - a}) ^ {(p + 1)/2}$ が答え

計算量は, どのステップも繰り返し二乗法を使って $\Theta(\log p)$ で計算できるので,

$\Theta(\log p)$ で平方根を求めることができました

参考文献

Cipolla のアルゴリズム – 37zigenのHP

平方剰余と基本的な問題 | 高校数学の美しい物語

ルジャンドル記号とオイラーの規準 | 高校数学の美しい物語

Cipolla's algorithm - Wikipedia

2021-09-14

Meissel–Lehmer algorithm (N以下の素数の数)

競技プログラミング

はじめに

Meissel–Lehmer algorithm ってな～～～に？

Library-Checker の問題のひとつ, Counting Primes

$N$ 以下の素数の個数を求めてください

$1 \le N \le 10^{11}$

これどうやってやるんですか？

ってことで調べたら出てきました

概要

思考をまとめるためにほとんどwiki写しただけです

素数計数関数を知っていますか？私は知っています

正の実数 $x$ について $\pi(x) = x$ 以下の素数の個数

という関数です

今回の問題は $\pi(N)$ を求めよという問題ですね

このままだと明らかにきついのでちょっと弱めたものを用意します

$p _ i = i$ 番目の素数として

$\pi(x, a) = ( [1, x ] \cap \mathbb{N} ) \setminus \displaystyle\bigcup _ {j = 1}^a p _ j \mathbb{Z}$

とします

？？？？となるかもしれませんが, エラトステネスの篩で $a$ 個目の素数まで篩ったときの素数候補って感じです

つまり, $1$ 以上 $x$ 以下の整数のうち $a$ 番目までの素数の倍数を除いたものです

さらに,

$\phi(x, a) = \left | \pi(x, a) \right |$

$P _ k(x, a) = \pi(x, a)$ のうち素因数をちょうど $k$ 個持っているものの数

とします

ここで, $k$ 個というのは $k$ 種類ということではありません例えば $12 = 2 \times 2 \times 3$ なので $3$ 個です

また, $P _ 0(x, a) = 1$ です ( $1$ は素因数を一つも持たない)

まず, 定義より

$\phi(x, a) = \displaystyle\sum _ {k = 0}^\infty P _ k(x, a)$

が成立することはわかると思います

更に,

$\pi(x) = P _ 1(x, a) + a$

が任意の $a$ について成り立ちます

$p _ a$ 以下の素数の倍数を全て篩ったあとに残ってるなかで素因数の数が $1$ のもの $= p _ a$ より大きい素数

よって $x$ 以下の素数の数 $= p _ a$ より大きい素数の数 $+ a$ となります

この二つの式を組み合わせると

$\pi(x) = a + \phi(x, a) - 1 - \displaystyle\sum _ {k = 2}^\infty P _ k(x, a)$

となります

ということで右辺に出てきたよくわからないのを二つ求められれば終わりです

というのが Meissel–Lehmer algorithm です

後は頑張ってください疲れました

もうひと頑張り

後は簡単ですね！とはいきません

なんか無限級数とか出されても困りますということでまだ続くんじゃ

まずは $\phi(x, a)$ の方から

漸化式のようなものを考えます小さいほうから篩っているので,

$\phi(x, a) = \phi(x, a-1) -$ ( $p _ a$ を最小の素因数とするような整数の数 )

というように考えられます

$p _ a$ を最小の素因数とするような整数の数をどうするかですが,

$p _ a$ の倍数だけを考えてみると, それらを $p _ a$ で割ったもののうち, $a - 1$ までの素数で篩ったものとなります (日本語にすると難しいですね)

ということでこれは $\phi \left( \dfrac{x}{p _ a} , a-1 \right)$ と一致します

したがって

$\phi(x, a) = \phi(x, a-1) - \phi \left(\dfrac{x}{p _ a}, a - 1 \right)$

という漸化式が得られます

次に $\displaystyle\sum _ {k = 2}^\infty P _ k(x, a)$ です

無限に足していくなんてやってられませんが, よく考えなくても $k$ には上限があります

もっと強くいきましょう

とある $a$ が存在して

$\displaystyle\sum _ {k = 2}^\infty P _ k(x, a) = P _ 2(x, a)$

となってくれたらとてもうれしいですね

素因数を三つ以上持つような合成数を全て篩い落とすには, $a$ をいくつにすればいいでしょうか

ということで $y = \sqrt[3]{x}$ として $a = \pi(y)$ としてみるといいんじゃないでしょうか

そんなわけで $P _ 2(x, a)$ が計算できればいいということになりました

定義からゴリゴリ行きましょう

$P _ 2(x, a) = \Big| \, n \, \big| \, n = p _ b p _ c, \, a \lt b \le c \, \Big |$

と書けますね

このとき $p _ b \le \sqrt{x}$ かつ $p _ c \le \dfrac{x}{p _ b}$ となるはずなので

$P _ 2(x, a) = \displaystyle\sum _ {b = a+1}^{\pi( \sqrt{x} )} \Big| \, n \, \big| \, n = p _ b p _ c, \, b \le c \le \pi \left(\dfrac{x}{p _ b} \right) \, \Big |$

と書けますこれは $c$ の数を数えればいいだけなので

$P _ 2(x, a) = \displaystyle\sum _ {b = a+1}^{\pi( \sqrt{x} )} \left( \pi \left(\dfrac{x}{p _ b} \right) - (b-1) \right)$

と書けて, 等差数列になってる項を括弧から出すと

$P _ 2(x, a) = \dbinom{a}{2} - \dbinom{\pi( \sqrt{x} )}{2} + \displaystyle\sum _ {b = a+1}^{\pi( \sqrt{x} )} \pi \left(\dfrac{x}{p _ b} \right)$

と書けました

はい, ようやく終わりですなんかまた $\pi$ が出てきてしまいましたがもともとの問題よりは小さくなっているので大丈夫です

まとめると

$\pi(x) = a + \phi(x, a) - 1 - P _ 2(x, a) \quad (y = \sqrt[3]{x}, a = \pi(y))$

$\phi(x, a) = \phi(x, a-1) - \phi \left(\dfrac{x}{p _ a}, a - 1 \right)$

$P _ 2(x, a) = \dbinom{a}{2} - \dbinom{\pi( \sqrt{x} )}{2} + \displaystyle\sum _ {b = a+1}^{\pi( \sqrt{x} )} \pi \left(\dfrac{x}{p _ b} \right)$

となりました

もうひと手間

いい感じに立式することはできましたが計算量はどうなっているでしょうか？

正直よくわからないですが, できればこんなに再帰させないでうまく計算してあげたいです

上で $y = \sqrt[3]{x}$ と書きましたが面倒なので $y$ は $\sqrt[3]{x} \le y \le \sqrt{x}$ である整数とします (最終的に $\sqrt[3]{x} \log ^ ? x$ とします)

まず $P _ 2(x, a)$ ですが, 素因数が二つしかないならなんかもっとうまくやれそうな気分になります

ここで登場するのは尺取り法です素因数が二つということでポインタを二つ使っていい感じにやります

まず, 前計算をします

エラトステネスの篩や線形篩を利用して $\dfrac{x}{y}$ 以下の素数のテーブルを作ります (計算量解析が面倒なので線形篩を使ったことにします)

ついでにそれぞれの整数に対して最小の素因数を求めておきます (これは後で使います)

さて, 定義より $P _ 2(x, a)$ として数えられるような数は, 素因数の小さいほうが $p _ a$ より大きいはずです

そのため素因数の大きいほうは $\dfrac{x}{p _ a}$ より小さいということになりますが, 前計算で $\dfrac{x}{y}$ 以下の素数は全て求めました

ということで, $p _ a \lt p \lt \dfrac{x}{p _ a}$ であるようなそれぞれの素数について, 素の素数がペアになるような素数の数を数えます

これは尺取り法を使えば素数の数に対して線形時間で数えることができますね

ということで $P _ 2(x, a)$ を求めることができました

さて, 残るは $\phi(x, a)$ の計算です, 再帰的に $\phi(m, b)$ を呼び出していくことになりますが

$m \le \dfrac{x}{y}$ のときは後でまとめて処理することにして今は考えないことにします

そうすると, 残りの呼び出されるものは $\phi \left(\dfrac{x}{k}, b\right) \quad k \lt y, b \le a$ です

これは $O(y \pi(y))$ 個ありますが, 素数定理より $\pi(x) \sim \dfrac{x}{\log x}$ という近似ができるので $O\left( \dfrac{y ^ 2}{\log y} \right)$ となります

ということでまぁいい感じに抑えられたということにしておきます

残っているのはさっきためておいた $m \le \dfrac{x}{y}$ の場合です

$\phi(m, b)$ は, 定義より $m$ 以下の自然数のうち最小の素因数が $p _ b$ より大きいものの数です

大きいものより小さいもののほうが数えやすそうですね

最小の素因数が $p _ b$ 以下であるものの数を数えることにします

ここで, だいぶ前に前計算しておいた最小の素因数を使いますうれしいことに, ちょうど $\dfrac{x}{y}$ 以下の整数について求めていました

ということで, $2$ 以上 $m$ 以下の整数の最小の素因数を列挙することはできます

$p _ b$ 以下であるものの数を数えますが, 求めたい $\phi(m, b)$ はたくさんあるので毎回愚直に数えるのではなく効率的に処理したいです

では, $m$ の昇順に見ていくことにしましょう

配列 $d$ を用意して, $d _ i = m$ 以下の整数のうち最小の素因数が $p _ i$ であるようなものの数とすると,

$\phi(m, b) = m - \displaystyle\sum _ {i = 1} ^ {b} d _ i$

と書けます区間和です

区間和なら累積和でできる？いいえ, 更新があるので Binary Indexed Tree (Fenwick Tree) を使いましょう

ということで, $\phi(m, b)$ は $m$ の昇順に見ていくことで合計 $\dfrac{x}{y} \log \dfrac{x}{y}$ で求めることができました

$\phi(x, a)$ 全体ではどうだったかというと $O\left( \dfrac{y ^ 2}{\log y} + \dfrac{x}{y} \log \dfrac{x}{y} \right)$ ですねよくわかりません

しかし, ここで $y = x ^ {1/3} \log ^ {2/3} x$ としてあげます

そうすると

$O\left( \dfrac{y ^ 2}{\log y} + \dfrac{x}{y} \log \dfrac{x}{y} \right) \\ = O\left( \dfrac{y ^ 2}{\log x} + \dfrac{x}{y} \log x \right) \\ = O\left( \dfrac{x ^ {2/3} \log ^ {4/3} x}{\log x} + \dfrac{x}{x ^ {1/3} \log ^ {2/3} x} \log x \right) \\ = O\left( x ^ {2/3} \log ^ {1/3} x \right)$

となります

ということで, 冒頭の

$N$ 以下の素数の個数を求めてください

というクエリについては計算量オーダー $O\left( N ^ {2/3} \log ^ {1/3} N \right)$ で解くことができました！

おわりに

マジで疲れました

英語記事しか見つからなかったので勉強しながら書きましたが, 理解が浅いので嘘が混じっているかもしれません優しく教えてください

まだ実装してませんが, 実装できたら追記しようと思います

ありがとうございました

追記: 通りました

f:id:suu_0313:20210915235416p:plain — https://judge.yosupo.jp/submission/60242

定数倍がよろしくない気がしますが, 実装の参考にどうぞ

参考文献

Meissel–Lehmer algorithm - Wikipedia

Counting primes in $$$\tilde{{\cal O}}(n^{2/3}\,\,)$$$ - Codeforces

2021-08-31

彩色数を求めるアルゴリズム

このスライドを読めばだいたいわかります

個人的に詰まったところメモ

最初の言いかえ

グラフ $G = (V, E)$ が $k$ 彩色可能 $\Leftrightarrow V$ を $k$ 個の独立集合で被覆できる

について, $k$ 個の独立集合は直交していなければいけないというような条件はない（それはそうで, $k$ 彩色で $k$ 色全部使わなくてもいいので同値）

独立に $k$ 個選ばないといけないと思っていたので, その後の

$f(S) = I(S) ^ k$

が重複しまくりやんけとなっていたが問題なかった

そのためこの $g(S)$ は彩色多項式 $P(G, k)$ とは異なる (ですよね？)

$P(G, k)$ を求めたかったらサイズの総和を $|V|$ で固定してあげればよくて, そうすると $f(S)$ はどうするんでしょうか

$I(S,n) =$ $S$ の部分集合で大きさ $n$ の独立集合の数を求められれば二次元DPにできそうな気分です

それはともかく後は頑張って $I(S)$ を求めればいい

独立集合の数え上げだが, $v \in S$ として

$v$ を要素として含むもの
$v$ を要素として含まないもの

のどちらかなので,

$I(S) = I(S - {v} - adj(v)) + I(S - {v})$

と書ける

ということで, あとはbitDPを書くだけですね

おしまい

2021-06-08

たまに見かける三種類の文字の二項演算について

AtCoder 競技プログラミング

はじめに

三種類の文字 $A, B, C$ があって二つの文字 $x, y$ について

$x = y$ のとき $x ( = y)$ を返す
$x \neq y$ のとき $x$ とも $y$ とも異なる文字を返す

というような二項演算が出てくる問題をたまに見かけます

こういう問題の解説を見ると~~よくわからない~~天才的な変換が書いてあることがしばしばあると思います。
それで納得できるならそれでいいのですが、そうもいかないという人に向けて~~全く参考になる気がしない~~自分がしている考え方を紹介しようと思います。（あくまでも「お気持ち」です）

お気持ち

この二項演算を $f$ とします

$f(A, A) = A, \quad f(A, B) = C$ となるような演算ということです。

さて、ここで登場する文字は三つだけなので、考えやすそうな数字に変換しましょう。

どんな順番でもいいですが、とりあえず

$A = 0$
$B = 1$
$C = 2$

とします。（問題が $R, G, B$ とかだったらそれに置き換えて考えてください）

こうすると、ここで求められる演算というのは

$f(0, 0) = 0$
$f(0, 1) = 2$
$f(0, 2) = 1$
$f(1, 0) = 2$
$f(1, 1) = 1$
$f(1, 2) = 0$
$f(2, 0) = 1$
$f(2, 1) = 0$
$f(2, 2) = 2$

こうなってほしいわけです。

さて、一気に考えるのは大変なのでとりあえずどちらも同じ数字のやつだけ考えてみます。

$f(0, 0) = 0$
$f(1, 1) = 1$
$f(2, 2) = 2$

突然ですが、二項演算と聞いてみなさん何を思い浮かべるでしょうか

だいたい最初に出てくるのは足し算・かけ算なのではと思っています。

ということでとりあえず足してみましょう。

とはいえそれだと二倍になってしまうので、足して 2 で割ってみましょう。

つまり、

$f(x, y) = \cfrac{x + y}{2}$

としてみます。こうするととりあえず今回考えている三つはみたしますね。

ついでに $f(0, 2) = 1$ なんかもみたしています。

しかし、このままでは $f(0, 1)$ とかで困ります。どうしましょうか

ここで、文字が三種類しかないということを生かしてみます

どうせ三つしかないので $\pmod 3$ の世界で考えてみるということです。

そうすると $2$ と $3$ は互いに素なので逆元が存在するはずで、

$2 \cdot 2 \equiv 4 \equiv 1 \pmod 3$

ということで

$2^{-1} \equiv 2 \pmod 3$

となります

よって

$f(x, y) \equiv 2(x + y)$

となり、

$f(0, 1) \equiv 2(0 + 1) \equiv 2 \pmod 3$

となってなんかいい感じです

そして、これはみたしてほしい関係をすべてみたしてくれます。

$2 \equiv -1$ なので $f(x, y) \equiv -x - y$ とも書けます。解説スライドでこの形の式を見かけた気がします。

と、こんな感じで二項演算の式を導くことができました。

割と自然な変換しかしていないつもりなのですが、どうでしょうか

例題

と、いうことでこの式で典型90 - 047 を考えてみます。

今回数えたいのは

${}^\forall i$ で $f(s_i, t_{i+k}) = c$ となるような $k$ の数です（ $c$ は定数）

さっきの式にばらして

$f(s _ i, t _ {i+k}) \equiv 2s _ i + 2t _ {i+k} \equiv c$

とします

ここで両辺に $2$ をかけると

$2(2s _ i + 2t _ {i+k} ) \equiv 2c$

$s _ i + t _ {i+k} \equiv 2c$

となって、 $t$ を移項して

$s _ i \equiv 2c - t _ {i+k}$

としてみます

ここで $- t _ {i+k}$ とはなんぞやということを考えてみると、

今三種類の文字を $0, 1, 2 \pmod 3$ として考えているのでそれをマイナスにしてみると

$-0, -1, -2 \equiv 0, 2, 1$

となりますね

E8さんの解説に登場した G→B、B→G という置き換えはこの式からも分かります。

また、 $c$ というのは何らかの定数ということにしているので

$u _ {i+k} \equiv 2c - t _ {i+k}$

と置いてみると

$s_i \equiv u _ {i+k}$

となって、ここまでくれば Z-algorithm などで共通する長さを求めてほげほげみたいなことをすればいいということがわかります。*1

おわりに

ということで、大した内容ではないですが天才じゃない人の思考過程にも多少は価値があるんじゃないかな～ってお気持ちで書いてみました。

いかがでしたか？

*1:ここから先はE8さんの解説などを見てください

2021-03-19

【作ってみた】二次元 Disjoint Sparse Table【なにこれ】

C++ 競技プログラミング

導入

訳あってSparse Tableについて調べていたのだが、こんな記事を見つけてしまった

この記事では「二次元のRmQクエリを高速に処理をする方法」として書かれているが、演算部分を抽象化することにより一般的な冪等半群のクエリを処理することができる二次元 Sparse Tableとなる。

二次元 Sparse Tableができるなら、二次元 Disjoint Sparse Tableもできるんじゃないか？

ひとまず検索

こういうのはだいたい先人が開発しているものである。別に車輪の再発名をすることに対して抵抗はないが、よりよい実装があるなら見てみたい。とりあえずGoogle先生に聞いてみることとしよう。

「二次元 Disjoint Sparse Table」や「2D Disjoint Sparse Table」など、それっぽい単語を入れてみる

出ない

なぜだ？

というかそもそも、2D Sparse Tableの時点でこどふぉのブログ一つ以外にそれっぽい記事が見当たらない。まぁ確かにこれが必要になるような問題ならみんな二次元セグ木を生やすのかもしれない...

それとも、そもそも二次元に拡張することが困難なのか？そんなことはないように思えるのだが...

書いてみる

仕方ないから1から自分で書くことにする。

上記のブログを参考にまずは2D Sparse Tableを書いてみて感覚をつかむ。これはお手本があるのですんなりと行った。tableの値をマージしていくのが面白かった。

ここからが本題、Disjointの方だ。tableの大きさはさっきと同じで $\log n \times n \times \log m \times m$ としてみる。

すると、 $\rm{table _ {0,i}}$ に $i$ 行目の要素についての 1D Disjoint Sparse Table を生やせばよさそうだ

$i$ にのループを回して生やしてみるここまで $O(nm\log m)$

そしたら次は生やしたDisjoint Sparse Tableを上手くマージして二次元にする

うーん、これはどうすればいいのだろう？ Sparse Tableに比べてマージの仕方が少しだけ複雑なので、脳内だけで考えると結構混乱してしまった。この辺は数学力が足りていないのだろう。

1D Disjoint Sparse Tableで区間長1の部分には引数に渡した元の配列の要素を入れるが、ここにさっき生やしまくったDisjoint Sparse Tableを使えばいいことに気づいた。これは行ける気がする。

とりあえずそれっぽいものが書けた。最後にクエリの処理をする方法を考えよう。

やはり縦横をそれぞれ二つに分けて合計四つの領域に分けるのがいいのだろう。しかし、Disjoint Sparse Tableは長さ1の区間は別処理をしなければいけない。

いい書き方が思いつかなかったので、丁寧に場合分けを書く。（といってもせいぜい四通りだ）

tableのどの部分を使うかについては、1D Disjoint Sparse Tableと全く同じ方法を使う。xorを取って、最上位bitの位置を求めるための前計算をする。

その前計算を保存しておく配列は、次元は増えても一次元ずつ処理していくので結局は変わらないわけだ。

そんなわけで完成してしまった。

template<typename T>
struct DisjointSparseTable2D{
  using F = function<T(T, T)>;
  const F f;
  vector<vector<vector<vector<T>>>> table;
  vector<int> lookup;

  DisjointSparseTable2D(const vector<vector<T>> &v, const F &f): f(f){
    const int n = v.size(), m = v[0].size();
    const int nb = 32 - __builtin_clz(n), mb = 32 - __builtin_clz(m);
    table.assign(nb, vector(n, vector(mb, vector<T>(m))));
    for(int rj = 0; rj < n; rj++){

      for(int cj = 0; cj < m; cj++) table[0][rj][0][cj] = v[rj][cj];

      for(int ci = 1; ci < mb; ci++){
        int shift = 1 << ci;
        for(int cj = 0; cj < m; cj += shift << 1){
          int t = min(cj+shift, m);
          table[0][rj][ci][t-1] = v[rj][t-1];
          for(int k = t-2; k >= cj; k--){
            table[0][rj][ci][k] = f(v[rj][k], table[0][rj][ci][k+1]);
          }
          if(m <= t) break;
          table[0][rj][ci][t] = v[rj][t];
          int r = min(t+shift, m);
          for(int k = t+1; k < r; k++){
            table[0][rj][ci][k] = f(table[0][rj][ci][k-1], v[rj][k]);
          }
        }
      }
    }

    for(int ri = 1; ri < nb; ri++){
      int rshift = 1 << ri;
      for(int rj = 0; rj < n; rj += rshift << 1){
        int rt = min(rj+rshift, n);
        for(int ci = 0; ci < mb; ci++){
          for(int cj = 0; cj < m; cj++){
            table[ri][rt-1][ci][cj] = table[0][rt-1][ci][cj];
            for(int k = rt-2; k >= rj; k--){
              table[ri][k][ci][cj] = f(table[0][k][ci][cj], table[ri][k+1][ci][cj]);
            }
          }
        }
        if(n <= rt) break;
        for(int ci = 0; ci < mb; ci++){
          for(int cj = 0; cj < m; cj++){
            table[ri][rt][ci][cj] = table[0][rt][ci][cj];
            int r = min(rt+rshift, n);
            for(int k = rt+1; k < r; k++){
              table[ri][k][ci][cj] = f(table[ri][k-1][ci][cj], table[0][k][ci][cj]);
            }
          }
        }
      }
    }
    const int b = max(nb, mb);
    lookup.assign(1 << b, 0);
    for(int i = 2; i < (1 << b); i++){
      lookup.at(i) = lookup.at(i>>1) + 1;
    }
  }

  T query(int x1, int x2, int y1, int y2) const {
    x2--; y2--;
    if(x1 >= x2){
      if(y1 >= y2) return table[0][x1][0][y1];
      int p = lookup[y1 ^ y2];
      return f(table[0][x1][p][y1], table[0][x1][p][y2]);
    }
    int np = lookup[x1 ^ x2];
    if(y1 >= y2) return f(table[np][x1][0][y1], table[np][x2][0][y1]);
    int mp = lookup[y1 ^ y2];
    T R1 = f(table[np][x1][mp][y1], table[np][x1][mp][y2]);
    T R2 = f(table[np][x2][mp][y1], table[np][x2][mp][y2]);
    return f(R1, R2);
  }
};

元の要素になる二次元配列と、マージする関数を引数に入れる設計です。 ei1333さんのライブラリの影響を多分に受けています（~~というかパクり~~）

使ってみる

さすがにバグが怖い

元となっている資料があるわけでもないし、何かよくわからない四次元配列を使った計算をしているんだから当然だ。

verifyができる問題を探すが、静的な二次元の区間minを要求される問題にそもそも出会ったことがない。Library-Checkerさんにもなさそう、困った。

見つけた～

さすがAOJ、何でもある

submit...

f:id:suu_0313:20210319011136p:plain — https://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=5299736#2

MLE！？

悲しい...

けど Case#1 は通っているので多分バグってはないと思います

追記:
テストケースをダウンロードして手元環境で動かしてみました。
出力結果はあっていたので大丈夫だと思います。

僕の実装を置いておきます

おわりに

ここまでありがとうございました。

この実装のバグや、もっといい実装、いい感じにverifyができる問題など見つけた方がいらっしゃいましたらご一報いただけると嬉しいです。

追記:
思考過程書きまくってたらろくに解説を書いていないことに気づきました。
そのうち加筆します。

2021-02-21

青コーダーになりました(現実)

初めに

色変記事的な言いたいことはだいたいこっちに書いたのでそちらを見てください

本文

SOMPO HD プログラミングコンテスト2021(AtCoder Beginner Contest 192) にて、ついに青コーダーになりました！！！！！！！！！！！！

f:id:suu_0313:20210220233848p:plain — Ratingグラフ

とても嬉しいです。正直上記の記事を書いてからたいして成長していないように思いますが、何とかここまでたどり着くことができました。

ここまで関わってきた皆様のおかげだと思います。本当にありがとうございます！！！

ギリギリ青色になれたとはいえ、当然まだまだ実力は足りていないと思います。これからも精進を続けて、もっともっと成長していきたいです。

とりあえずの目標としては、ABC以外でも結果を出せるようにすることでしょうか

f:id:suu_0313:20210220235709p:plain — https://rating-history.herokuapp.com/rating.html?handle=Suu0313&number=3

このように、ABC以外ではひどい結果となってしまっているので...

また、春休みに入ったのでゆきこやこどふぉなどにも手を出してみたいなと思っています。

おわりに

青コーダーになるまで競プロを続けてこれたのは、いつもアドバイスをくれたり絡んでくれたこれを読んでくれている皆さま、また、もちろんAtCoderさまのおかげです。

ほんとうにありがとうございます！これからもよろしくお願いします。