エラトステネスの篩 よくない実装と計算時間

エラトステネスの篩

エラトステネスの篩は、素数を列挙するアルゴリズムです。詳しい説明は省きますが、10までの素数を列挙したいなら、trueで初期化したサイズ10の配列を作成し、2の倍数をfalseに、3の倍数をfalseに、、、という操作を繰り返すことで、最後までtrueで残った整数が素数とわかるというものです。

計算量

(エラトステネスの篩には様々な高速化のTIPSがありますが、)エラトステネスの篩をまともに実装した場合の漸近的な計算量は $O(n \log \log n)$ です。しかし、実装を少し間違えると計算量が $O(n \log n)$ になってしまったりします。

合成数でも篩を走査してしまう場合

(2021/08/22修正)

import math
def eratosthenes_O_nlogn(n):
  sieve = [True] * (n + 1)
  for i in range(int(math.sqrt(n)) + 1):
    if i < 2:
      sieve[i] = False
    else:
      for j in range(2, n//i + 1):
        sieve[i * j] = False
  return sieve

このコードは、サイズ $n+1$ のリストを篩として用意することで、インデックス $i$ へのアクセスで、 $i$ が素数であるかどうかを判定できるようにしています。for文は $n$ の平方根まで回せば十分です($i^2>n$ のときに break という処理で実現した方がいいとは思いますが、今回は平方根まで回すということを明示するためにmath.sqrtを使用しています)。 $i<2$ のもの、つまり $i = 0, 1$ については素数ではないものとして、 False を代入し、次のループへ進みます。その他の$i \ge 2$ については、その倍数の篩に、Falseを代入していきます。

～計算量の解析～

さて、このコードでも十分高速に素数の列挙ができるのですが、漸近的な計算量は $O(n \log n)$ です。それは、合成数であってもその倍数の篩にFalseを代入する処理を行ってしまっているためです。篩の走査でFalseが代入される回数を簡単に解析してみます。とはいっても立式は簡単で、整数 $i$ について、$n/i$ 回Falseが代入されるので、 $i$ を $[2, \sqrt{n}]$ の範囲で足してやればいいだけです。

$$
\frac{n}{2} + \frac{n}{3} + \frac{n}{4} + \frac{n}{5} + \cdots + \frac{n}{\sqrt{n}} = \sum_{i=2}^{\sqrt{n}} \frac{n}{i} = n \sum_{i=2}^{\sqrt{n}} \frac{1}{i}
$$

$n$ は $i$ に関係ないので、Σの外へ出すことができます。総和の部分は調和級数の形になっており、逆数の和が、 $y=\frac{1}{x}$ の積分が $\log x$ であることなどから、$\log n$ で漸近的に抑えられることが知られています。今回は、$\sqrt n$ までしか足さないので、結局、$O(n \log \sqrt n)$ で抑えられることになります。ここで、$\log$ の性質や $O$ の性質を使うと、$\log$ の中の $\sqrt{}$ は消すことができるため、$O(n \log \sqrt n) = O(n \log n)$です。

$$ n \sum_{i=2}^{\sqrt{n}} \frac{1}{i} = O(n \log \sqrt n) = O(\frac{1}{2} n \log n) = O(n \log n)
$$

以上から、合成数でも篩を走査してしまうコードだと、計算量が $O(n \log n)$ になってしまうことがわかりました。

まともに実装したエラトステネスの篩

さて、先ほどのコードでは合成数でも篩を走査してしまっていました。しかし、そのせいで、計算量が $O(n \log \log n)$ ではなく、 $O(n \log n)$ になってしまうということを解析しました。合成数では篩を走査しないようにすれば、その問題は回避できます。

(2021/08/22修正)

import math
def eratosthenes(n):
  sieve = [True] * (n + 1)
  for i in range(int(math.sqrt(n)) + 1):
    if i < 2:
      sieve[i] = False
    elif sieve[i]:
      for j in range(2, n//i + 1):
        sieve[i * j] = False
  return sieve

新たなコードでは、篩を走査する条件にsieve[i]を追加して、 $i$ がすでに素数でないことが判明している場合は篩を走査せず、次のループへ進むようになっています。

～計算量の解析～

篩の走査でFalseの代入が実行される回数を簡単に解析してみます。先ほどと同様、整数 $i$ について、$n/i$ 回Falseが代入されるのですが、今回の $i$ は素数のみとなります。つまり、素数の逆数の和を考えることになります。

$$
\frac{n}{2} + \frac{n}{3} + \frac{n}{5} + \frac{n}{7} + \cdots = n \left(\frac{1}{2} + \frac{1}{3} + \frac{1}{5} + \frac{1}{7} + \cdots \right) = n \sum_{p \le \sqrt{n}}\frac{1}{p}
$$

$n$までの素数の逆数の和が $\log \log n$ で抑えられることは過去の偉人によって示されているそうなので( https://ja.wikipedia.org/wiki/メルテンスの定理 )、オーダーの計算ができます。今回のコードは$\sqrt{n}$までループを回すことに注意します。

$$
n \left(\frac{1}{2} + \frac{1}{3} + \frac{1}{5} + \frac{1}{7} + \cdots \right) \
= n \sum_{p \le \sqrt{n}}\frac{1}{p} = O(n \log \log \sqrt n) = O(n (\log \log n + \log \frac{1}{2})) = O(n \log \log n) \
$$

以上から、まともに実装したエラトステネスの篩の計算量は $O(n \log \log n)$ となります。

おまけ ループをnの平方根で切らない場合

$n$ の平方根でループを切らない場合を考えてみます。ただし、まともに実装したとき同様、合成数では篩を走査しません。

(2021/08/22修正)

def eratosthenes_O_nloglogn(n):
  sieve = [True] * (n + 1)
  for i in range(n):
    if i < 2:
      sieve[i] = False
    elif sieve[i]:
      for j in range(2, n//i + 1):
        sieve[i * j] = False
  return sieve

このコードの計算量は、$O(n \log \log n)$ で、まともに実装したエラトステネスの篩と同じ計算量になります。

～計算量の解析～

今までの解析では、実は、ループを平方根で打ち切ることによって発生する $O(\sqrt{n})$ の項は、最終的には消えていました。そのため、ループを平方根で打ち切らなくても、漸近的な計算量は変わらないことになり、$O(n \log \log n)$ です。

計測

さて、上では3つのエラトステネスの篩のようなアルゴリズムを考えました。合成数でも篩を走査するもの、まともに実装したもの、平方根でループを打ち切らないものの3つで、計算量はそれぞれ $O(n \log n)$、$O(n \log \log n)$、$O(n \log \log n)$ です。

それぞれのアルゴリズムについて、$n=10^1$ ～ $10^6$からほぼ等間隔に50個の値をとり、それぞれの $n$ について100回ずつ試行した実行時間の平均をグラフにしてみました。実行環境はGoogle Colaboratoryです。

(2021/08/22修正)

縦軸は実行時間で単位は秒、横軸は $n$ です。

やはり $O(n \log n)$ の実装は $O(n \log \log n)$ のどちらの実装よりも漸近的には遅そうだという結果になりました。また、計測した範囲では、ループを平方根で打ち切らない場合は打ち切る場合よりも2倍程度時間がかかりそうということがわかりました。

エラトステネスの篩を実装するとき、ループを平方根で打ち切るのを忘れるよりも、合成数でも篩を走査してしまうことのほうがよくないということが、解析からも計測結果からも示されていそうです。

付録

計測に用いたコード

import numpy as np
from tqdm import notebook
import timeit

number = 100

x = np.linspace(10, 10**6, num=50, dtype=int)

eratosthenes_time = []
eratosthenes_O_nloglogn_time = []
eratosthenes_O_nlogn_time = []
for i in notebook.tqdm(x):
  eratosthenes_time.append(timeit.timeit(lambda: eratosthenes(i), number=number)/number)
  eratosthenes_O_nloglogn_time.append(timeit.timeit(lambda: eratosthenes_O_nloglogn(i), number=number)/number)
  eratosthenes_O_nlogn_time.append(timeit.timeit(lambda: eratosthenes_O_nlogn(i), number=number)/number)

(eratosthenes_time, eratosthenes_O_nloglogn_time, eratosthenes_O_nlogn_time)

グラフの描画に用いたコード

import numpy as np
import matplotlib.pyplot as plt

plt.plot(x, eratosthenes_time, label='origin(nloglogn)')
plt.plot(x, eratosthenes_O_nloglogn_time, label='nloglogn')
plt.plot(x, eratosthenes_O_nlogn_time, label='nlogn')
plt.legend()