エラトステネスの篩とMeissel-Lehmer Algorithm

エラトステネスの篩

競プロのために、いろいろなアルゴリズムを自力で実装してストックすることを始めた。 今のところ、以下をライブラリとして Git に push している。

• Next Permutation
• Union Find Tree
• Binary Search：二分探索
• Cumulative Sum：累積和
• Eratosthenes Seive：エラトステネスの篩

今日はエラトステネスの篩を実装していたが、思ったように速度が出ない。Python の遅さも相まって、コンテストで実用に耐えうるのはせいぜい$N=10^7$といったところである。（手元の環境でも AtCoder のコードテストでも 700ms 程度かかる）

素数の全列挙が必要なケースでは、エラトステネスの篩で戦わざるを得ないが、単純に素数の個数を求めるだけであれば、より優れたアルゴリズムがあるようなので、調べていた。これがかなり難しく、なかなか理解が追いつかないので、思考過程も含めてここにアウトプットしておく。

Meissel-Lehmer Algorithm

なぜ高速化を志したかというと、Library-Checker での入力が$N=10^{11}$で、話にならんレベルで高速化が必要だからである。

「こんなん C++専用すぎて Python じゃ通せないだろ」と思って AC 例をあさっていたら、なんとちゃんと AC している人がいる。 強い。これは強い。アルゴリズムの可能性をひしひしと感じた。

幸いにも使用しているアルゴリズムを引用元を記載のうえコードを提出している神のような Submission があったので、勢いいさんで調べてみたもののゴリゴリの整数論でノックアウト寸前である。

それが今回の Meissel-Lehmer Algorithm ということなのだが、ちょっと書ききれるかわからない。というか理解できるかもわからないのが現状である。

現時点でわかっていること

まず、前提となるエラトステネスの篩は、対象となる整数列を小さい方から順番に見ていって、2 の倍数、3 の倍数、5 の倍数……という調子に数を消していくアルゴリズムである。 順番に消していけば、$n$の倍数を消すステップに遷移した際に、必ず$n$は素数である。なぜなら、ある素数$p$を考えた際に、$p$以下でかつ$p$と互いに素な数は$p-1$までのすべての整数であるからだ。それまでにふるい落とされていないということは、$n$は$n-1$までのすべての数と互いに素である。

ここで、ある自然数$x$に対して$x$以下かつ$x$と互いに素な自然数の個数を$\varphi(x)$とする。これをオイラーのファイ関数あるいはトーシェント関数などと一般的に呼ぶ。

$p$を素数とすると、$\varphi(p^k)=p^k-p^{k-1}=p^k(1-\frac{1}{p})$という性質や、互いに素な数$m,n$に対して、$\varphi(mn)=\varphi(m)\varphi(n)$という性質を持つ。

前者は、$p^k$の素因数が$p$のみであることを考えれば、$p^k$の素因数でないものは$p$の倍数に限定される。それは当然、$p^k$を$p$で割れば求まるため、導出できる。

後者は、乗法性と呼ばれるがこれが少し難しい。まず、$m\times{n}$の表を作って考える。

$m/n$	$1$	$2$	$3$	$\cdots$	$t$	$\cdots$	$m$
$1$	$1$	$2$	$3$	$\cdots$	$t$	$\cdots$	$m$
$2$	$m+1$	$m+2$	$m+3$	$\cdots$	$m+t$	$\cdots$	$2m$
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
$n-1$	$(n-2)m+1$	$(n-2)m+2$	$(n-2)m+3$	$\cdots$	$(n-2)m+t$	$\cdots$	$(n-1)m$
$n$	$(n-1)m+1$	$(n-1)m+2$	$(n-1)m+3$	$\cdots$	$(n-1)m+t$	$\cdots$	$nm$

$t$と$m$と互いに素な数だとする。 $\varphi$関数の定義から、1 行あたり$t$となりうる数は$\varphi(n)$個ある。このとき、$t$で構成される列は$n$の完全剰余系であり、そこに含まれる$n$と互いに素である数の個数は$\varphi(n)$となる。

したがって、各$t_1,t_2,\cdots,t_{\varphi(m)}$に対して、$\varphi(n)$個ずつ$m,n$と互いに素となる数があるので乗法性が成り立っている。

完全剰余系における互いに素である数の個数については要補足。今後の予定。

Prime counting function

$x$以下に含まれる素数の個数を返す関数を素数計数関数、または$\pi$関数という。 $\pi(10)=4$のように表す。$10$以下の素数は、$2,3,5,7$の$4$つである。

Partial seive function

ここで$\pi(x)$は、既往のエラトステネスの篩を用いれば求めることが可能であるが、$\pi(x)$をもっと計算負荷の低い何か別の関数への操作で得ることができれば、計算量を落とすことが可能である。

そこで、部分ふるい関数（partial seive function）という概念を導入する。 これは、$\displaystyle\phi(x,a):=\left|([1,x]\cap\mathbb{N})\backslash\bigcup_{j=1}^{a}p_j\mathbb{Z}\right|$で定義される集合である。

式だけ書くと意味不明なので、言葉で解釈すると閉区間$[1,x]$にある自然数を、$a$番目までの素数でふるった残りの数である。

ここで、$a$番目の素数とは、$\{p_1,p_2,p_3\cdots\}=\{2,3,5\cdots\}$というように$2$から小さい順に素数にインデックスを割り付けたときの番号を指す。

$k$th partial seive function

さらに、これを計算するために$k$次の部分ふるい関数を定義する。

これは、$P_k(x,a)$で表し、$\phi(x,a)$のうち、素因数をちょうど$k$個だけ含む数の集合を表す。

$k$次部分ふるい関数の具体例

例として、$\phi(10,1)$を考えてみる。

$10$以下の素数は$2,3,5,7$であり、$a=1$であるので、$1$番目の素数$2$でふるった残りとなることが分かる。

したがって、以下となる。

さらに、このうち素因数を$1$つも含まないものを部分ふるい関数で表すと$P_0(10,1)=\{1,\cancel3,\cancel5,\cancel7,\cancel9\}=1$となる。

また、素因数をただ$1$つ含むもの、すなわち素数それ自身は、$P_1(10,1)=\{\cancel1,3,5,7,\cancel9\}=3$である。

部分ふるい関数と$k$次部分ふるい関数の関係性

ここで$k$次のふるい関数で得られる集合は、次数が異なればそれぞれの集合は排他的であるので、以下が言える。

また、以下の式も成り立つ。これは、定義を考えれば自明であろう。

さらに、上の$(1)$式を$P_1$について解けば以下も成り立つことがわかる。

したがって、$(2)$および$(3)$を連立して以下が言える。

力尽きた…

ちょっと重くなってきたので、いったん記事を切ることにする。 そう簡単には「かんぜんにりかいした」段階にすら至らない重さだ。 終わりが見えない…