素数判定プログラム -- 3

さて，下請け関数をコンパイラ言語で書いたらどうなるかということで，私が一番とっつきやすい FORTRAN（雀百まで踊り忘れず）を使って，以下のようなプログラムを書いた。

function is_prime(x)
    implicit none
    integer(8)  mx, x, i
    logical is_prime
    if (x == 2 .or. x == 3) then
            is_prime = .true.
      return
    else if (x <= 1 .or. mod(x, 2) == 0 .or. mod(x, 3) == 0) then
      is_prime = .false.
      return
    end if
    mx = int(sqrt(float(x)))
    do i = 5, mx, 6
      if (mod(x, i) == 0) then
        is_prime = (x == i)
        return
      else if (mod(x, i + 2) == 0) then
        is_prime = (x == i + 2)
        return
      end if
    end do
    is_prime = .true.
end function is_prime

function next_prime(x)
    implicit none
    integer(8) x, i, next_prime
    logical is_prime
    x = x + 1
    next_prime = 0
    do i = x, 9223372036854775807_8 ! 巨大整定数の宣言 2^63 - 1
        if (is_prime(i)) then
            next_prime = i
            return
        end if
    end do
end function next_prime

これを，

$ f2py -m prime_lib -c is_prime.f90

で

prime_lib.cpython-38-darwin.so

というシェアードライブラリを作る。

それを利用するには，まずインポートして

import prime_lib

prime_lib.next_prime() のようにする。

from time import time
start = time()
k = 9223372036854775782
print(f"next_prime({k}) = {prime_lib.next_prime(k)}")
print(time() - start)

実行結果は

next_prime(9223372036854775782) = 9223372036854775783
0.00039386 秒で答えが出る（小数点以下3桁の秒数なんてほぼゴミなので無視して良いレベル）。

Python で書いたプログラムでは 197.028 秒だったから，50万倍くらいの速度だ。

Julia で書いたプログラムでは 9.343 秒だったから，2万4千倍くらいの速度だ。

Julia は大健闘しているものの，やはりコンパイラ言語には太刀打ちできない。

素数判定プログラム -- 2

さて，前の記事で Python による素数判定プログラムを評価したが，今回は Python プログラムを Julia に移植した結果を書く。

まず，2 を特別にして，あとは奇数で試し割りをするプログラムである。

Julia は色々な変数タイプが使えるが，符号付き 64　ビット整数の Int64 を使うことにする。

使用法は，x::Int64 のように変数のアノテーションを行う（初出の一回だけでよいようだ）。

他の型からの変更は Int64() のように行う。

function is_prime0(x::Int64)
    if x == 2
        return true
    elseif x <= 1 || x % 2 == 0
        return false
    end
    mx::Int64 = Int64(floor(sqrt(x)))
    for i::Int64 in 3:2:mx
        if x % i == 0
            return false
        end
    end
    return true
end

2 と 3 を別扱いにして 6n ± 1 で試し割りを行うバージョンが以下のプログラム。

function is_prime(x::Int64)
    if x == 2 || x == 3
        return true
    elseif x <= 1 || x % 2 == 0 || x % 3 == 0
        return false
    end
    mx::Int64 = Int64(floor(sqrt(x)))
    for i::Int64 in 5:6:mx
        if x % i == 0
            return x == i
        elseif x % (i + 2) == 0
            return x == i + 2
        end
    end
    return true
end

これを利用して，与えた数より大きい最小の素数を返す関数を書く。

function next_prime(X::Int64)
    X += 1
    answer::Int64 = 0
    for i::Int64 in X:9223372036854775783
        if is_prime(i)
            answer = i
            break
        end
    end
    return answer
end

実行時間の測定は 対象とする文の前に '@time' を付けるだけでよい。

@time println(next_prime(Int64(9223372036854775782))) # 9223372036854775783, 9.343510s

Python の場合

next_prime(9223372036854775782) は 197.028s だったので，ほぼ 20 倍速いということになった。

 

素数判定プログラム -- 1

原作者が最初に書いた関数。
本人も認めるとおり，速度的には超マズイ。

def is_prime000(x: int) -> bool:
    if x <= 1:
        return False
    for i in range(2, x):
        if x % i == 0:
            return False
    return True

試し割りはその数の平方根までで十分。
そこまでの数で割りきれなければ，素数である。

import math
def is_prime00(x: int) -> bool:
    if x <= 1:
        return False
    for i in range(2, (math.floor(math.sqrt(x)) + 1)):
        if x % i == 0:
            return False
    return True

原作者はここまで努力はしたものの，もっと速い方法（コンパイラを使う）を探ることにしてしまった。
しかし，よいアルゴリズムを使わないと，コンパイラでもカバーしきれないのはよく知られていることだ。もう少しプログラムを改良しよう。

素数は 2 だけが偶数で，それ以外は奇数。
よって，2 を特別扱いし，その他の素数候補の試し割りを行うのは奇数だけ。ということを使ったのが次のプログラム。

def is_prime0(x: int) -> bool:
    if x == 2:
        return True
    elif x <= 1 or x % 2 == 0:
        return False
    for i in range(3, int(x**0.5 + 1), 2):
        if x % i == 0:
            return False
    return True

さらに改良を加えよう。
5 以上の素数は 6n ± 1 の形をしている。
従って，2，3 を特別扱いし，その他の素数候補の試し割りを行うのは  6n ± 1 だけ。

def is_prime(x: int) -> bool:
    if x in [2, 3]:
        return True
    elif x <= 1 or x % 2 == 0 or x % 3 == 0:
        return False
    for i in range(5, int(x**0.5 + 7), 6):
        if x % i == 0:
            return x == i
        elif x % (i + 2) == 0:
            return x == i + 2
    return True

さて，4 つのプログラムの速度を評価しよう。
1000000007 は素数であるが，is_prime000(), is_prime00(), is_prime0(), is_prime() の計算時間を比較すると，80 秒，0.003 秒，0.002 秒，0.001 秒程度である。しかし，is_prime000() 以外は，誤差の範囲内でしかない。

from time import time
start = time(); print(is_prime000(1000000007), time() - start) # True 80.07880806922913
start = time(); print(is_prime00(1000000007), time() - start) # True 0.0031898021697998047
start = time(); print(is_prime0(1000000007), time() - start) # True 0.001580953598022461
start = time(); print(is_prime(1000000007), time() - start) # True 0.001107931137084961

しかし，92233720368547841 が素数かどうかを判定するのは，is_prime() でも 15 秒かかる。

start = time(); print(is_prime(92233720368547841), time() - start) # True 15.08924412727356

is_prime0() が 20 秒，is_prime00() が 40 秒ほどである。

start = time(); print(is_prime0(92233720368547841), time() - start)  # True 20.365103244781494
start = time(); print(is_prime00(92233720368547841), time() - start) # True 40.88630700111389

上の 1000000007 の素数判定にかかる計算時間との関係から，もし is_prime000() で判定するとしたら 80:0.003 = x:40 つまり, x = 80/0.003 * 40  = 1066666.667 秒 = 12.3 日 かかる。

さて，素数判定プログラムができれば，一つの利用法は，与えられた数より大きい最小の素数を求めるプログラムを書くこと。

def next_prime(x):
    x += 1
    answer = 0
    for i in range(x, 9223372036854775783 + 1):
        if is_prime(i):
            answer = i
            break
    return answer

start = time(); print(next_prime(10000000000000000), time() - start) # 10000000000000061, 7.773s

start = time(); print(next_prime(92233720368547757), time() - start) # 92233720368547841, 15.040s

大きな数を素数判定しようとすると，かなりの時間が掛かることになる。

start = time(); print(next_prime(9223372036854775782), time() - start) # 9223372036854775783, 197.028s

ここではじめてコンパイラ言語の使用を検討することにするが，その前に，Julia でプログラムしてみよう。

次の 13 日の金曜日はいつか？

ツェラーの公式を利用して列挙する

def Zeller(y, m, d):
    w = ['Sat', 'Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri']
    if m == 1 or m ==2:
        m += 12
        y -= 1
    C = y // 100
    Y = y % 100
    return w[(d + 26*(m+1)//10 + Y + Y//4 -2*C + C // 4) % 7]

for y in range(2021, 2030):
    for m in range(1, 13):
        if Zeller(y, m, 13) == 'Fri':
            print(y, m)

2021 8
2022 5
2023 1
2023 10
2024 9
2024 12
2025 6
2026 2
2026 3
2026 11
2027 8
2028 10
2029 4
2029 7