算数とプログラミングで遊ぼうの会¶

素数判定第三回: 素数判定の高速化¶

前回は合成数判定を高速化しました. 今回は素数判定そのものを高速化します.

補足: どのくらい速度が改善できるのか?¶

先に書いておきましょう. 前回のプログラムは整数$N$の大きさに比例する計算量があります. 特に$N$に比例する計算時間が必要です. 今回紹介する簡単な改善では計算時間は$\sqrt{N}$に比例するレベルに落ちます. これがどのくらい劇的な改善か, これまた具体的に計算して確認します.

次のコード片で10**nは$10^n$です.

In [1]:
import math
for n in [1,10**2,10**3,10**4,10**5,10**10,10**20]:
    print()
    print(f"n = {n}")
    print(f"sqrt({n}) = {math.sqrt(n)}")
    print(f"sqrt(n)/n = {math.sqrt(n)/n}")

n = 1
sqrt(1) = 1.0
sqrt(n)/n = 1.0

n = 100
sqrt(100) = 10.0
sqrt(n)/n = 0.1

n = 1000
sqrt(1000) = 31.622776601683793
sqrt(n)/n = 0.03162277660168379

n = 10000
sqrt(10000) = 100.0
sqrt(n)/n = 0.01

n = 100000
sqrt(100000) = 316.22776601683796
sqrt(n)/n = 0.00316227766016838

n = 10000000000
sqrt(10000000000) = 100000.0
sqrt(n)/n = 1e-05

n = 100000000000000000000
sqrt(100000000000000000000) = 10000000000.0
sqrt(n)/n = 1e-10

リスト内の各nに対して$\sqrt{n}$と$\frac{\sqrt{n}}{n}$を計算しています. この$\frac{\sqrt{n}}{n}$はまさに計算にかかる時間がどのくらい圧縮されるかを表しています. 実際に計算するとわかるように, 数が小さいうちは大したことはありません. 大きくなればなるほど計算量(所要時間)が$\sqrt{n}$になる影響が大きくなります.

もしあなたがプログラムを組んだ経験がないなら, こう言われただけでは実感が湧かないでしょう. これから実際に素数判定でどのくらい計算時間に影響が出るのか確認します.

結論の関数¶

解説はあとにしてまずは計算量(計算時間)が大きく減らせる関数を作り, 前回の関数と計算時間を比較します. 特に計算開始・終了に対するCPU時間の差を表示します.

前回の関数¶

In [2]:
def is_prime_ver2(n):
    b = True
    i = 2
    while i < n:
        if n%i == 0:
            b = False
            break
        i = i+1
    return b

高速化関数¶

In [3]:
def is_prime_ver3(n):
    b = True
    i = 2
    while i*i <= n:
        if n%i == 0:
            b = False
            break
        i = i+1
    return b

どこが違うかというとwhileの終了条件で, 前者はi < n, 後者はi*i <= nです.

高速化版の挙動確認¶

まずは正しく素数判定できるか確認しましょう.

In [4]:
for n in range(2, 40):
    if is_prime_ver2(n) == is_prime_ver3(n):
        print(f"{n}: The same result!")
    else:
        print(f"{n}: The wrong result!")

2: The same result!
3: The same result!
4: The same result!
5: The same result!
6: The same result!
7: The same result!
8: The same result!
9: The same result!
10: The same result!
11: The same result!
12: The same result!
13: The same result!
14: The same result!
15: The same result!
16: The same result!
17: The same result!
18: The same result!
19: The same result!
20: The same result!
21: The same result!
22: The same result!
23: The same result!
24: The same result!
25: The same result!
26: The same result!
27: The same result!
28: The same result!
29: The same result!
30: The same result!
31: The same result!
32: The same result!
33: The same result!
34: The same result!
35: The same result!
36: The same result!
37: The same result!
38: The same result!
39: The same result!

39以下の整数に対して既存の関数と結果が一致しました. 最低限の確認ができたので先に進みます.

補足: 簡潔なプログラムの紹介¶

実はリスト内包表記とPython組み込みのall()関数を使うとプログラムが簡潔に書けます. このコンテンツはプログラミング初心者向けなので, よけいな混乱を招かないようにここでは詳しく説明しません. こうしたプログラムを書くためのキーワードとして, 関数型言語・関数プログラミングを紹介するに留めます.

In [5]:
[is_prime_ver2(n) == is_prime_ver3(n) for n in range(2,40)]
Out[5]:
[True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True,
 True]

これにall()を作用させます.

In [6]:
all([is_prime_ver2(n) == is_prime_ver3(n) for n in range(2,40)])
Out[6]:
True

大きな数で計算時間を比較する¶

実際に計算するとすぐわかるように, ある程度大きな数にならないと計算時間の差を体感できません. 一方で数が大きすぎると計算が終わらなくなります. 体感できるいい塩梅を探るのも兼ねていくつか計算してみましょう.

n=100000から10倍して確認¶

順に確認します. 素数判定の高速化を確認したいのでまずは適当に素数を見つけましょう. 適当に次のような計算をすれば候補は簡単に見つかります.

In [7]:
for n in range(100000, 100100):
    if is_prime_ver3(n):
        print(n)

100003
100019
100043
100049
100057
100069

ここではn=100003を採用します.

n=100003: 十万¶

In [8]:
%%time
n = 100003
is_prime_ver2(n)

CPU times: total: 46.9 ms
Wall time: 39 ms
Out[8]:
True

ver3版も調べます.

In [9]:
%%time
is_prime_ver3(n)

CPU times: total: 0 ns
Wall time: 0 ns
Out[9]:
True

結果への注意¶

2022-02時点での私の手元の環境ではis_prime_ver3()は早く終わり過ぎてうまく計算できていません. 特にis_prime_ver2()は31.2ms程度である一方, is_prime_ver3()は0nsと出ました. この時点で高速化版が実際に$10^6$倍以上速くなっています.

コンピューターに計算させればすぐ終わる¶

いまnは十万のオーダーです. もちろん人間ではこんなに一瞬で計算できません. この圧倒的な計算速度がコンピューターに計算させる醍醐味です.

10n: 百万¶

is_prime_ver2()の計算はmsかかっていて既にかなり重いです. あなたはもしかしたら「このくらい体感だと一瞬なのに?」と思ったかもしれません. しかし私はこのコンテンツを作る試行錯誤の中で何十回・何百回とプログラムを実行しています. 特に高速化版も走らせているため, msしかかからない処理でも重く感じるのです.

そこで前回の低速版の測定はここで打ち止めにします. (私の手元の環境で)is_prime_ver3()の測定結果がまともに出てくるところまで数を持ち上げます. オーダーを一つ上げるとn = 1000003が素数です.

低速版¶

In [10]:
%%time
is_prime_ver2(1000003)

CPU times: total: 312 ms
Wall time: 310 ms
Out[10]:
True

高速版¶

In [11]:
%%time
is_prime_ver3(1000003)

CPU times: total: 0 ns
Wall time: 996 µs
Out[11]:
True

結果への注意¶

相変わらず高速化版はすぐに計算が終わります. まだ一秒かかっていないので低速版も手計算よりは比べるのも失礼なほど速いです. 低速版はnが10倍になったことで計算時間も大まかに10倍になった点にも注意してください.

100n: 千万¶

さらにオーダーを一つ上げます. ここからは高速版だけ調べます. もちろんあなた自身がいろいろ試すときには低速版でも試してみてください.

今度の素数はn = 10000019です.

In [12]:
%%time
is_prime_ver3(10000019)

CPU times: total: 0 ns
Wall time: 1.99 ms
Out[12]:
True

結果への注意¶

ぜひipynbで低速版を実際に実行してみてください. 私の実行環境ではここまでとは違って明らかに待たされる感覚があります. 一方is_prime_ver3()はいまだに差が取れないほど速いです.

1000n: 一億¶

さらに数を10倍しました. 今度の素数はn = 100000007です.

In [13]:
%%time
is_prime_ver3(100000007)

CPU times: total: 15.6 ms
Wall time: 5 ms
Out[13]:
True

結果への注意¶

私の手元ではとうとう1.99msといったミリ秒の時間がかかるようになりました.

10000n: 十億¶

さらに数を10倍しました. 今度の素数はn = 1000000007です.

In [14]:
%%time
is_prime_ver3(1000000007)

CPU times: total: 15.6 ms
Wall time: 14 ms
Out[14]:
True

結果への注意¶

私の計算環境ではそろそろ確実にmsオーダーの時間がかかりはじめました. もっと増やしてみましょう.

100000n: 百億¶

さらに数を10倍しました.

In [15]:
%%time
is_prime_ver3(10000000019)

CPU times: total: 62.5 ms
Wall time: 52.1 ms
Out[15]:
True

結果への注意¶

私の計算環境でも明確にmsオーダーです. 理屈の上ではnが10倍になっても所要時間は3倍, 多くても5倍程度です. この揺れはいろいろな事情に依存する部分で, 実際には何度も実行して平均や分散を調べる統計処理が必要です.

もっと試してもいいのですが, きりがないのでこのコンテンツとしてはいったんここで打ち切ります. ぜひ興味に応じてあなた自身で試せるだけ・試したいだけ試してみてください.

計算結果のまとめ¶

では計算結果をまとめます. 前回のis_prime_ver2()は100003の時点で既にmsオーダーでした. しかしis_prime_ver3()はその十万倍でようやくmsオーダーの時間です.

理論的な計算時間のオーダー解析¶

実は次のようになっています.

  • 前回作った関数is_prime_ver2(): nが10倍になると計算時間も10倍になる.
  • 今回作った関数is_prime_ver3():
    • is_prime_ver2()よりも速い.
    • nが10倍になると計算時間は3(=$\sqrt{10}$)倍程度にしかならない.

このnに関する挙動の違いはもちろんプログラムの組み方によります. まずはプログラムの挙動が何に由来するかを考え, そのあとにプログラムの改善に対する数学的観点を説明しましょう.

プログラムの再掲¶

前回の関数¶

def is_prime_ver2(n):
    b = True
    i = 2
    while i < n:
        if n%i == 0:
            b = False
            break
        i = i+1
    return b

高速化関数¶

def is_prime_ver3(n):
    b = True
    i = 2
    while i*i <= n:
        if n%i == 0:
            b = False
            break
        i = i+1
    return b

違いはwhileのi < nとi*i <= nだけ¶

本当にこれだけです. これで上記のように猛烈な速度差が出ます.

慣れていないとこのコードの意図は読み取れないかもしれません. 気分としては次のような処理をしています.

import math
sqrt_n = int(math.sqrt(n))
while i <= sqrt_n:

つまりループを回す回数が$n$から$\sqrt{n}$に激減しています. これが, まさにこれだけが速度差の原因で, オーダーが1上がるごとに新関数での計算時間が3倍程度にしか増えない理由です. せっかくなので$\sqrt{10}$をプログラムで計算して値を確認します.

In [16]:
import math
print(math.sqrt(10))

3.1622776601683795

いちいちプログラムで計算しなくても, $\sqrt{9} = 3$よりは大きく, $\sqrt{16} = 4$よりは小さいという暗算レベルの数的感覚も重要です.

なぜ$\sqrt{n}$まででいいのか?¶

もちろん数学的な理由があります. ある合成数$N$が因数$a$を持つとしましょう. 当然$b = N/a$も整数で$1$より大きいです.

ここで$a < b$を仮定すると, $ab = N$だから$a^2 < ab = N$で$a < \sqrt{N}$が得られます. 一方$b^2 > ab = N$によって$b > \sqrt{N}$も得られます. つまりある因数が得られればそれと対になる$\sqrt{N}$より大きい値も自動的に得られます. 特にいまは素因数分解したいわけではなく, 素数か合成数か判定したいだけなので, $\sqrt{N}$より小さい範囲で因数を調べれば十分なのです.

i*i <= nの意図¶

ここにはプログラミング上の工夫がたくさんあります. 少し長くなりますが私が思いつく理由をいくつか挙げておきます. プログラミングに慣れていないと何を言っているかわからないと思うので, 「そんな事情があるのか」と流し読みしても構いません. むしろこれを軸にプログラミングを勉強してもらうのも一手です.

整数と浮動小数点の型の問題¶

数学的にはwhile i <= math.sqrt(n)で構いませんし, Pythonで計算するときにも問題にはなりません.

しかし他のプログラミング言語, 特にコンパイルが必要な言語では問題になりえます.

何が問題かというと値の型です. ふつうiは整数の中で回します. 一方, 一般にsqrt(n)は数学的に実数で, プログラミング上は浮動小数点数です. 整数と浮動小数点数を不等式で比較しようとすると, 値(または変数)の型が合わないと言われてエラーになるプログラミング言語があります. つまり, そもそもとしてwhile i <= math.sqrt(n)と書くとプログラムが実行できないのです. 整数同士で比較するためにわざと実数の中での不等式i <= math.sqrt(n)と等価な整数の中での不等式i*i <= nを使っています.

都度計算の問題・演算量の問題¶

上のプログラムではwhile i <= math.sqrt(n)と書かず, あえて次のように書きました.

import math
sqrt_n = int(math.sqrt(n))
while i <= sqrt_n:

ここでは$\sqrt{n}$を先に計算してその結果を使ったのがポイントです.

プログラミング言語や最適化の状況にもよるのですが, 実はwhile i <= math.sqrt(n):と書くと, whileの条件判定をするごとにいちいち平方根を取る処理が走ります. ループの変数iは都度変わる一方, 平方根の$\sqrt{n}$は一回計算したらもう変わらないため, 何度も計算するのを避けるべく上のような書き方にしています.

平方根を取る処理の重さも問題です. 人間の手間で考えてもわかるように, 平方根を取る処理は整数のかけ算に比べて非常に重たい処理です. つまり都度のi <= math.sqrt(n)よりも, i*i <= nの方が遥かに軽いのです.

演算量の問題¶

上の注意にはさらなる問題があります. もしnが非常に大きい場合, 最初に$\sqrt{n}$を計算した上で切り上げ処理などで整数化しておき, i <= sqrt_nのように書いた方が速い可能性があります.

条件判定i*i <= nでは都度i*iを計算しています. 整数同士の積だから軽めの処理だとは言え, 何億回もの計算の果てでは塵も積もれば山になります. どちらの方が速いかは実際に計算して実測しないとわかりません. 状況によっては違いは統計的な誤差の範囲でおさまる可能性があります.

「推測するな. 計測せよ.」¶

上で実際に計算して実測しないとわからないと書きました. これはプログラムを書く上で非常に重要です. 節タイトルの標語としてよく知られています. 先程も書いたようにプログラミング言語にも依存する部分もあれば, 同じプログラミング言語でもコンパイル時の最適化オプション設定にも依存する場合があります.

何をどういう設定で考えればいいかは状況次第なのです. 状況に合わせてきちんと考え, 実測して何がいいかを判定する必要があります.

補足: いろいろな高速化の方法¶

速度単純にプログラムを書き直すにもいろいろな対処法があります. それこそ「遅いけれども書き直すコストが莫大だから諦める」という手もあります. 書き直すにしても「Pythonではなくそもそもが高速なC言語で書き直す」という手があります. 何をどう書くかにも依存するので一概にどのくらい速くなるとは言えませんが, 全く同じロジックで書いたプログラムが一万倍以上速くなることもあります

  • 【Python】for文の速度を上げよう～C言語で高速化～.

多少なりともプログラムを書ける人はたいてい一つだけではなく, 複数のプログラミング言語でプログラムを書けます. もちろん自然言語と同じで, 複数のプログラミング言語でプログラムが書けると言っても特別に得意な言語もあれば, 簡単な読み書きしかできない言語もあります. はじめのうちはPythonだけで構いませんが, 慣れてきたらぜひ他の高速なプログラミング言語にアタックしてみてください. 単に言語を変えるだけで爆速になる可能性はあります.

微分方程式を解くとき, 特に物理や工学で研究用にプログラムを書くときはよくC・C++・Fortranが使われます. 最近はJuliaの採用も増えているようですが, こと数値計算への応用には非常に強い三言語です. 現状スパコン上で動かすことを考えるとこの三言語しかほぼ選択肢がないと聞いたこともあります.

細かい話はともかく, 2-3言語くらいは眺めてみると見える景色も広がります.

補足: 平方根でどのくらい値が変わるか?¶

数学系プログラミングの勉強を進めていくと, 必ず計算量のオーダーという概念に出くわします. その中でデータ量や処理量nに応じて, 必要な計算の量がO(n^k)なのかO(2^n)なのかが問題になります. 後者は指数オーダーと呼ばれ, 絶対に避けるべき計算です. 前者はkが小さければ小さいほど嬉しくなります. ここではkが小さいとどのくらい嬉しいか, 数値の肌感覚を持ってもらうために具体的に計算して値を確認します.

特に$n$と$\sqrt{n}$の違いをいろいろな$n$で計算します. 見やすくするために平方根は整数化した上でスペースでパティングしています.

In [17]:
import math
for n in [10**k for k in range(0, 10)]:
    print(f"n      = {str(n).rjust(10, ' ')}")
    print(f"sqrt n = {str(int(math.sqrt(n))).rjust(10, ' ')}")
    print()

n      =          1
sqrt n =          1

n      =         10
sqrt n =          3

n      =        100
sqrt n =         10

n      =       1000
sqrt n =         31

n      =      10000
sqrt n =        100

n      =     100000
sqrt n =        316

n      =    1000000
sqrt n =       1000

n      =   10000000
sqrt n =       3162

n      =  100000000
sqrt n =      10000

n      = 1000000000
sqrt n =      31622

特に最後のn = 1000000000を比較してください. 確認する数に30000倍もの違いがあります. これがnを固定したときの二つの関数の速度の違いです. 30000倍のレベルで速度差があるため, 前者が1秒かかっているときに, 後者が雑な処理では)検出不可能なレベルの短時間になったのです. ただルートを取っただけでこれだけの差が出る凄みを味わってください.

今回の高速化のまとめ¶

今回も数学的な事情を使って素数判定に必要な数を大幅に減らして高速化しました.

問題に応じていろいろな高速化法があります. 今回, 数学的には定義を詳しく検討して定義の条件を緩める形で数学的に厳密に対応しました. 実際にはPythonという原理的に処理が遅いプログラミング言語ではなく, C言語のように処理が速い言語を使う手もあります.

応用からするともっと凄い高速化法もあります. その名もずばり, 確率的素数と呼ばれる概念, ミラーラビン法による確率的素数判定と呼ばれる手法・アルゴリズムがあります. そしてミラーラビン法自体, 実は合同式やフェルマーの小定理のような厳密な数学を背景にしています. よいアルゴリズムを作るプログラミングの問題解決のためにそれ相応の数学が必要な局面があり, 素数判定はまさにそうした問題なのです.

次回予告¶

プログラミングへの誘いとしての素数判定コンテンツはいったんここで終わります. 次回はさらなるプログラミング学習のための情報・コンテンツを紹介します. あなたの興味・関心にしたがってどんどん勉強してください.