唯物是真 @Scaled_Wurm

プログラミング(主にPython2.7)とか機械学習とか

素因数分解とかエラトステネスの篩(ふるい)とかのメモ

素数を求めたり素因数分解するのは競技プログラミングでたまに出てきます
計算量とか詳細をあまり知らなかったので基本的なアルゴリズムについて調べてみました

アルゴリズムや数学についてはあまり詳しくないので間違いがあったら指摘してください
ランダウの記号\(O(\cdot)\)とか使ってますが理解はゆるふわです

まずは基本的なところから始めていきます
以下では正の整数についてだけ考えます

素数

\(1\)とそれ自身でしか割ることができない整数(\(1\)は含まれない)
$$ 2, 3, 5, 7, 11, 13, 17, 19, \dots $$

素因数分解

整数を素因数(約数になる素数)の積で表す
$$ 60 = 2^2 \times 3 \times 5 $$

試し割り(\(\sqrt n\)以下のすべての整数)

ある数を素因数分解する為の簡単なアルゴリズム
いろいろな整数で割ってみて、割り切れたらその素数で何回割れるかをカウントしていく

ある数を\(n\)とした時に\(\sqrt n\)以上の素因数は最大でも\(1\)個しか含まれていないので\(\sqrt n\)以下のすべての整数で割れるか確かめれば十分

例えば\(n=22\)なら\(\sqrt 20 \approx 4.69\)なので\(2, 3, 4\)で割り切れるか確かめればよい
\(2\)で\(1\)回割りきれて残りの数は\(11\)
\(11\)は\(3, 4\)で割り切れない
なので\(22\)の素因数分解の結果は\(22 = 2 \times 11\)となる

\(\sqrt n\)以下のすべての整数について割るので\(O(\sqrt n)\)
素因数の数は\(2\)の冪乗の時が最大なので、素数をカウントする部分は\(O(\log n)\)

なので(時間)計算量はだいたい\(O(\sqrt n)\)

\(\sqrt n\)を計算するときに誤差のおそれがあるので実際には\(\sqrt n + 1\)まででループしたほうがよいと思います

import math
import collections
def trial_division_sqrt(n):
    prime_count = collections.Counter()
    
    for i in xrange(2, int(math.sqrt(n)) + 2):
        while n % i == 0:
            n /= i
            prime_count[i] += 1
    if n > 1:
        prime_count[n] += 1
    
    return prime_count

エラトステネスの篩(ふるい)

エラトステネスの篩は素数を見つけるためのアルゴリズムとしてよく知られているものです
複数の素数を列挙する場合には、試し割りによる方法よりも高速に処理できます

どのようなアルゴリズムかという以下の画像がわかりやすいです
このアルゴリズムは\(n\)以下のすべての数について素数かどうかを表した表を作成します

まず最初にすべての数を素数だと仮定して初期化します
\(2\)から順に数を見ていって、もしその数が素数だったら\(n\)以下のその倍数を素数でないとマークしていきます
たとえば2が素数だとしたら\(4, 6, 8, \dots\)、3が素数だとしたら\(6, 9, 12, \dots\)を素数でないとしてマークします
この手順を\(n\)まで繰り返します

Sieve of Eratosthenes animation.gif
"Sieve of Eratosthenes animation". Licensed under CC BY-SA 3.0 via Wikimedia Commons.

次に計算量について考えます
たとえば\(2\)の倍数を確認するとき\(\frac{n}{2}\)回処理を行うので、\(i\)番目の素数\(p_i\)については\(\frac{n}{p_i}\)回処理します
すべての素数について考えると、\(n\)以下のすべての素数の個数を\(\pi(n)\)と置くと
$$\frac{n}{2} + \frac{n}{3} + \frac{n}{5} + \frac{n}{7} + \dots = \sum_i^{\pi(n)} \frac{n}{p_i}$$

素数の逆数の総和\(\sum_i^{\pi(n)} \frac{1}{p_i}\)は漸近的に\(\ln \ln n\)に近づくらしいので、\(n\)までの素数を列挙するときのエラトステネスの篩の計算量は以下のようになります
$$O(\sum_i^{\pi(n)} \frac{n}{p_i}) = O(n \sum_i^{\pi(n)} \frac{1}{p_i}) = O(n \log\log n)$$

偶数はすべて除いたり、ある素数\(p\)の倍数を確認するときには\(p^2\)以降だけをみればよいとかいろいろ高速化できるみたいです(\(p^2\)未満の\(p\)の倍数はすでにふるわれている)

ちなみに上で書いた試し割りの方法で\(n\)までの素数を求める計算量はたぶん\(O(n\sqrt n)\)ぐらいなので結構な差になります

import itertools
def eratosthenes_sieve(n):
    table = [0] * (n + 1)
    prime_list = []
    
    for i in xrange(2, n + 1):
        if table[i] == 0:
            prime_list.append(i)
            for j in xrange(i + i, n + 1, i):
                table[j] = 1
    
    return prime_list

試し割り(素数だけを使う)

上で書いた試し割りでは\(\sqrt n\)までのすべての整数で割りましたが、\(\sqrt n\)までの素数だけで割るともっと速くなります
試し割りというと主にこっちの方法を指すっぽいです(?)

\(n\)以下の素数の個数は\(\pi(n) = \frac{n}{\ln n}\)と近似できるらしいです

試し割りの時には\(\sqrt n\)までの素数がわかればよいので計算量は\(O(\frac{\sqrt n}{\log n})\)
あらかじめ\(\sqrt n\)までの素数を列挙しておくのにエラトステネスの篩を使うと、前処理の計算量として\(O(\sqrt n \log\log n)\)かかります

ちなみに\(n\)までの素数列挙を試し割り(素数)でするときの計算量よりも、エラトステネスの篩の方が速いらしいです
検索して出てきた結果によると、素数の試し割りで求める方法は\(\frac{n\sqrt n}{(\log n)^2}\)

performance - Sieve of Eratosthenes versus trial division time complexity - Stack Overflow
http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf

import math
import collections
def trial_division_prime(n, prime_list):
    prime_count = collections.Counter()
    
    for i in prime_list:
        if i * i > n:
            break
        while n % i == 0:
            n /= i
            prime_count[i] += 1
    if n > 1:
        prime_count[n] += 1
    
    return prime_count

エラトステネスの篩の結果を使った素因数分解

エラトステネスの篩と同様の処理を行って、表にある数が素数かどうかだけではなく、ある数の素因数を一つ記録しておきます(メモリに余裕があれば重複を除いたすべての素因数でもよい)
素因数分解するときには、表でその数が素数かどうかを見て、素数でない時には記録してある素因数の一つで割ってその値について表がどうなっているかを順番にたどっていきます
こうすれば含まれる素因数以外で割ることがなくなるので高速化できます

前準備の表を作る部分の計算量はエラトステネスの篩と同様に\(O(n \log\log n)\)
実際に素因数分解で表をたどる部分は、素因数の数は\(2\)の冪乗の時に最大なことから\(O(\log n)\)となります

import itertools
import collections
def prime_factor_table(n):
    table = [0] * (n + 1)
    
    for i in xrange(2, n + 1):
        if table[i] == 0:
            for j in xrange(i + i, n + 1, i):
                table[j] = i
    
    return table

def prime_factor(n, prime_factor_table):
    prime_count = collections.Counter()
    
    while prime_factor_table[n] != 0:
        prime_count[prime_factor_table[n]] += 1
        n /= prime_factor_table[n]
    prime_count[n] += 1
    
    return prime_count

実際に素因数分解の速度がどれぐらい違うか試してみた

\(n=10^5, 10^6, 10^7, 10^8\)としたときに\([n - 10000, n]\)の範囲の数を素因数分解する時間を比較しました
以下の時間の単位は秒です

\(n\) 試し割り(\(\sqrt n\)まで) 試し割り(\(\sqrt n\)までの素数) エラトステネスの篩の結果を使った素因数分解
\(10^5\) 0.22 0.05 0.04
\(10^6\) 0.65 0.08 0.27
\(10^7\) 2.04 0.14 2.76
\(10^8\) 6.36 0.28 32.41

結果を見ると一万個程度の素因数分解をするだけなら試し割り(エラトステネスの篩で列挙した素数だけ使う)が一番速かったです

次に\(n=10^5, 10^6, 10^7, 10^8\)としたときに\(n\)までのすべての整数を素因数分解する時間を比較しました

\(n\) 試し割り(\(\sqrt n\)まで) 試し割り(\(\sqrt n\)までの素数) エラトステネスの篩の結果を使った素因数分解
\(10^5\) 1.61 0.48 0.31
\(10^6\) 45.76 7.03 3.36
\(10^7\) 125.94 34.98

この場合はエラトステネスの篩の結果を使った素因数分解の方が速かったです

エラトステネスの篩の結果を使った素因数分解が、少ない数を素因数分解するときに遅いのは前処理のエラトステネスの篩が重いからです
この前処理の計算量は\(O(n \log\log n)\)ですが、試し割り(\(\sqrt n\)までの素数)の方の前処理の計算量は\(O(\sqrt n \log\log n)\)なので、あまり多くない数を素因数分解するときには前処理の時間が支配的になります

蛇足

素因数分解の結果やエラトステネスの篩的なアルゴリズムは約数の個数などを求めるのにも役立ちます

約数の個数と総和

ある数\(n\)が以下のように素因数分解できるとき
$$n = \prod_{i=1}^r {p_i}^{a_i}$$
約数の個数や総和は以下の約数の\(x\)乗の総和の式で求められます
$$\sigma_x(n) = \prod_{i=1}^r \sum_{j=0}^{a_i} {p_i}^{jx} = \prod_{i=1}^r (1 + {p_i}^x + {p_i}^{2x} + \dots + {p_i}^{a_i x})$$
\(x=0\)の時約数の個数
$$\sigma_0(n) = \prod_{i=1}^r (a_i + 1)$$
\(x=1\)の時約数の総和
$$\sigma_1(n) = \prod_{i=1}^r (1 + {p_i}^1 + {p_i}^{2} + \dots + {p_i}^{a_i})$$

def num_divisors(prime_count):
    num = 1
    
    for prime, count in prime_count.iteritems():
        num *= (count + 1)
    
    return num

def sum_divisors(prime_count):
    sum = 1
    
    for prime, count in prime_count.iteritems():
        temp = 1
        temp_sum = 1
        for i in xrange(count):
            temp *= prime
            temp_sum += temp
        sum *= temp_sum
    
    return sum

また素因数分解して上の公式を使う以外にも、エラトステネスの篩的な方法を使ってまとめて計算することもできます
\(n\)以下のすべての整数\(i\)についてループして、\(i\)の倍数について素数かどうかをマークするのではなく、約数の個数を求めるなら\(+1\)ずつ、約数の和を求めるなら\(+i\)ずつ加えてカウントしていけばよいです
各整数\(i\)について\(\frac{n}{i}\)回ずつ処理を行うので計算量は調和級数の和の\(n\)倍になるので次のようになります
$$O(\sum_{i=1}^n \frac{n}{i} = \frac{n}{1} + \frac{n}{2} + \frac{n}{3} + \frac{n}{4} + \frac{n}{5} + \dots + \frac{n}{n}) = O(n \log n)$$

def num_divisors_table(n):
    table = [0] * (n + 1)
    
    for i in xrange(1, n + 1):
        for j in xrange(i, n + 1, i):
            table[j] += 1
    
    return table

def sum_divisors_table(n):
    table = [0] * (n + 1)
    
    for i in xrange(1, n + 1):
        for j in xrange(i, n + 1, i):
            table[j] += i
    
    return table

オイラーのトーシェント関数

前に以下の記事でも書きましたが、\(n\)以下の共通の素因数を持たない(互いに素な)整数の個数を求めるオイラーのトーシェント関数というものがあります

TopCoder SRM 617 Div1 --- 1350->1268 ――オイラーのトーシェント関数 - 唯物是真 @Scaled_Wurm
オイラーのトーシェント関数は\(n\)の\(i\)番目の素因数を\(p_i\)としたとき次のようになります
$$\varphi(n) = n \prod_{i=1} \left(1-\frac{1}{p_i}\right)$$

def totient(prime_count):
    num = 1
    
    for prime, count in prime_count.iteritems():
        num *= (prime - 1)
        for i in xrange(1, count):
            num *= prime
    
    return num

これも同様にエラトステネスの篩的に計算ができます
素数かどうかの表とは別に、\([1, n]\)の数値そのものを入れた表を作って、素数\(p_i\)の倍数について\(\left(1-\frac{1}{p_i}\right)\)をかけていけば計算できます
計算量はエラトステネスの篩と同じく\(O(n \log\log n)\)

def totient_table(n):
    table = range(n + 1)
    
    for i in xrange(2, n + 1):
        if table[i] == i:
            for j in xrange(i, n + 1, i):
                table[j] *= float(i - 1) / i
    
    return map(lambda x: int(x + 0.01),  table)

参考

他にもエラトステネスの篩の高速化ができたり、アトキンの篩とかいろいろなアルゴリズムがあるみたいです

追記: 線形時間なエラトステネスの篩(2014-11-17)

\(O(n)\)でエラトステネスの篩(?)をやる方法がTwitterで紹介されていたのでメモっておく

以下のアルゴリズムでは、合成数の最小の素因数を順番に求めていく
現在注目している数を\(i\)としたとき、\(i\)の最小の素因数よりも小さな素数を\(p\)とおく
このとき\(p\times i\)の最小の素因数は\(p\)となる
各々の合成数はエラトステネスの篩の場合と違って、最小の素因数についてしか見られない
つまり最小の素因数のテーブルは合成数の個数回しか更新されない
なので計算量は\(O(n)\)になるはず

def linear_sieve(n):
    min_prime_factor = range(n + 1)
    prime = []
    for i in xrange(2, n + 1):
        if min_prime_factor[i] == i:
            prime.append(i)
        for p in prime:
            m = p * i
            if p <= min_prime_factor[i] and m <= n:
                min_prime_factor[m] = p
            else:
                break
    return prime

ちなみに単純な実装で実験したら普通のエラトステネスの篩の方が速くなりましたorz
……なんか間違えてるのかな?

100000000までの素数を見つけるのを10回繰り返した時の合計時間

linear: 565.265198167
eratosthenes: 319.974106872

計測に使ったコード

import timeit

linear_setup = """
def linear_sieve(n):
    min_prime_factor = range(n + 1)
    prime = []
    for i in xrange(2, n + 1):
        if min_prime_factor[i] == i:
            prime.append(i)
        for p in prime:
            m = p * i
            if p <= min_prime_factor[i] and m <= n:
                min_prime_factor[m] = p
            else:
                break
    return prime
"""
linear = 'linear_sieve(100000000)'

eratosthenes_setup = """
def eratosthenes_sieve(n):
    table = [0] * (n + 1)
    prime_list = []
    
    for i in xrange(2, n + 1):
        if table[i] == 0:
            prime_list.append(i)
            for j in xrange(i + i, n + 1, i):
                table[j] = 1
    
    return prime_list
"""
eratosthenes = 'eratosthenes_sieve(100000000)'

print 'linear:', timeit.timeit(linear, linear_setup, number=10)
print 'eratosthenes:', timeit.timeit(eratosthenes, eratosthenes_setup, number=10)