互いに素である a,b の組を効率良く列挙しよう

互いに素である 2 つの自然数を効率的に列挙する方法について、2 通りのアルゴリズムを紹介する。

Stern-Brocot tree を使い列挙する

この木がどんなものであるかは、Wikipedia を読むのが一番手っ取り早い。 とはいっても、英語版にはそのもののページがあるが、日本語版には無いようだ。しかし、ファレイ数列のページでこの木について触れられているので、そこのリンクを貼っておく。

• Stern–Brocot tree - Wikipedia, the free encyclopedia
• ファレイ数列 - Wikipedia # 隣接する分数と中間数

ユークリッドの互除法を逆算する要領で列挙する

最大公約数を求める際にユークリッドの互除法が使われるが、それを逆の手順で計算する事によっても、互いに素である組を網羅する事が出来る。つまり、

$${\rm gcd}(a,\: b) = \begin {cases} a & (b = 0) \\ {\rm gcd}(b,\: a \bmod b) & (otherwise) \end {cases} \\ {\rm gcd}(b,\: a \bmod b) = {\rm gcd}(b,\: a+bn)$$

である事を利用し、${\rm gcd}(1,\: 1) = 1$ から始めて逆算を行えば、全ての互いに素である a,b の組を列挙可能となる。

C++ と Scheme による実装

まずは C++ で書いたものから。

そして Scheme で書いたもの。

2 つのアルゴリズムの違い

まず 1 つ目に、再帰関数に渡さなければならない引数の数が違う。引数の数が違えば、再帰呼び出しにおけるスタックの消費量が違い、また push/pop の要不要によって処理速度も変わってくる。その為、Stern-Brocot tree によるものは、ユークリッドの互除法を逆の手順にしたものよりも、メモリを食い速度も若干遅くなる。

2 つ目に、特定の区間にある既約分数を列挙出来るかどうかといった差異がある。Stern-Brocot tree を用いたものでは、enum_coprimes_with_stern_brocot_tree(N, 0, 1, 1, 1) として、$[\frac01,\: \frac11]$ に存在する既約分数を列挙可能だが、ユークリッドの互除法を逆の手順にしたものでは、簡単には実現出来ない。因みに、Stern-Brocot tree をもちいたアルゴリズムで $[\frac13,\: \frac12]$ に存在する既約分数を列挙するには、enum_coprimes_with_stern_brocot_tree(N, 1, 3, 1, 2) とすれば良い。これは他の区間でも同様である。

ちょっとした応用

これまで、互いに素となる a,b の組、つまり ${\rm gcd} (a,b) = 1$ となる a,b の組を列挙をする方法を書いてきたが、少し応用すれば ${\rm gcd}(a,b) = n$ となる a,b の組を列挙する事も出来る。これは単純に、互いに素となる a,b の組を列挙し an,bn としても計算しても良いが、以下のようにする事でも可能だ。

こうしてやれば、互いに素である a,b の組を生成する都度 n を掛けるより、若干効率が良い。

Project Euler - Problem 432 を強引に素早く解こう

問題文はこちら

• http://projecteuler.net/problem=432
• http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20432

$S(n,m)$ を計算する上手い方法が分からないので、素直に $\varphi(n \cdot i)$ を全て求めて足し上げる事で答えを出す事とする。しかし、$\varphi(n \cdot i)$ のように n, i を掛けてから考えると面倒なので、下のように分けて考える。

$$\varphi(n \cdot i) = \varphi(n)f(i)$$

このようにしておけば、$S(n, m)$ の定義は

$$S(n,m) = \sum_{i=1}^m \varphi(n \cdot i) = \varphi(n) \sum_{i=1}^m f(i)$$

となり、$f(i)$ に集中して考える事が出来る。$f(i)$ の求め方は i を素因数分解して

$$i = \prod_{j=1}^r {p_j}^{k_j}$$

と与えられるならば、以下のように計算する事が出来る。

$$f(i) = \prod_{j=1}^r g(p_j^{k_j}) \\ g(p^k) = \begin{cases} \varphi(p^k) = (p-1)p^{k-1}, & (\mbox{if } gcd(p, n) = 1) \\ p^k, & (otherwise) \end{cases}$$

尚、ここで $\varphi(p^k)$ を求める際、速度的な理由から除算を使ってはいけない。

S(510510, 10^6) を求めよう

1 から N までの φ を適用した値を欲しい時、エラトステネスの篩とほぼ同様の手順で計算する事が出来る。f の場合も同様に可能であり、その考え方を用いて $S(510510, 10^6)$ を求めるプログラムを C++ で書くとこのようになる。

このプログラムを簡単に説明すると、$S(510510, N = 10^6)$ を計算する際、まず全ての $f(i)$ の値を保存するテーブル coefs[N+1] = { 0, 1, 1, ..., 1 } を生成し、 そして一旦 N までの素数列を全て列挙してから、その各素数の倍数について $f(i)$ を段階的に計算していく事で、$\varphi(510510)$ に掛けるべき係数を求めるものとなっている。因みに先述したように、count_coefs() における for 文の形から、エラトステネスの篩とほぼ同様の手順で $f(i)$ を計算している事が分かる。

$S(510510, N = 10^9)$ でも 50s で求められているのでそこそこ高速ではあるが、上記のプログラムでは N 個の uint32_t の要素を持つテーブルを生成しなければならないので、N = 10^9 ならおよそ 4GB 必要となる。その上 N = 10^11 となると $f(i)$ が uint32_t で収まらない可能性があるので、uint32_t でなく uint64_t でテーブルを持つ必要があり、その場合およそ 800GB 必要となってしまい、後述する工夫が必要となる。

S(510510, 10^11) を求めよう

先程のプログラムには大きな欠点が 2 つある。

1. N = 10^10, 10^11 では、篩で使うテーブルサイズが大きくなり過ぎて実行出来ない
2. N までの素数列を全て列挙しなければならない

1 つ目の欠点については、区間篩を使い、小さな区間に分けて $f(i)$ を列挙する事で解決出来る。区間篩を使う際、篩に掛けるテーブルサイズに、キャッシュメモリに乗るだけのサイズを指定すれば、DRAM へのアクセスが減り高速化を期待する事も出来るので、一石二鳥にもなる。

2 つ目の欠点を解消するには、$f(i)$ を列挙する篩で使う素数列を、$\sqrt{N}$ までのものに限定すれば良い。$\sqrt{N}$ までの素数列を用いて、篩に掛ける要領で $f(i)$ を部分的に計算した時、coefs[] の中に全くアクセスされず 1 のまま変動しないものが出てくるが、その要素を coefs[x] とすれば、その x こそが $\sqrt{N}$ より大きな素数を意味している。値が変動しない要素というのは、$\sqrt{N}$ までのどの素数でも割り切れない数を指しており、それは $\sqrt{N}$ より大きな素数に他ならないからだ。

そんなわけで、区間篩を使い、$\sqrt{N}$ までの素数列で篩うようにしたプログラムがこちら。区間篩を使うと、各区間を並列計算出来るようになるので、マルチスレッド化も序でに施してある。

このプログラムの補足をしておく。$\sqrt{N}$ までの素数列を small-primes とし、$\sqrt{N}$ より大きな素数列を big-primes と呼ぶとする。small-primes でそれぞれの区間を篩に掛け、big-primes を徐々に見つけていくのだが、この見つけた big-primes を用いて篩を掛ける必要は無い。何故なら big-primes における各素数 p の、N までの全ての p の倍数に対する $f(p \cdot n)$ の総和は、$\varphi(p) S(510510, \lfloor \frac{N}{p} \rfloor)$ と簡単に計算可能だからだ。また $\frac{N}{p} < \sqrt{N}$ と充分に小さく、$S(510510, 1), \: ..., \: S(510510, \lfloor \sqrt{N} \rfloor)$ をキャッシュしておけば、メモリに一度アクセスし $\varphi(p)$ を掛けるだけで計算出来る。

ぶん回して 4 分なのでそこそこ。以前にも同じような感じで区間篩とマルチスレッドを用いて、1e10 までの素数がいくつあるか、2 秒前後で数えるプログラムを書いた事があったので、問題文を読んですぐに、わりと短めな処理時間で解けると分かった。

というわけで

• 区間篩でメモリアクセスの局所化やっぱりおすすめ
• 区間篩 + OpenMP による並列処理もやっぱりおすすめ
• MathJax おすすめ

でした。