素数判定・素数生成プログラム
$1$ と自身以外に約数をもたない自然数を 素数 (prime number) といいます。言い換えると、$2$ 個の約数をもつ自然数は素数です。たとえば、$7$ の約数は $1$ と $7$ の $2$ 個だけなので素数です。$1$ の約数は $1$ だけなので素数ではありません。$2$ は $1$ と $2$ の $2$ 個の約数をもつので素数です。すなわち 2 は最小素数です。$3$ 以上の素数はすべて奇数なので、$2$ は唯一の偶数の素数 (偶素数) です。
素数判定のアルゴリズムはシンプルです。ある自然数 $N$ が与えられたとき、$N$ を $2$ から順に $N/2$ を超えない数まで割って、一度も割り切れなければ素数です。たとえば、$N=11$ のとき、$2,\ 3,\ 4,\ 5$ で割ってみると、一度も割り切れないので素数です ($N/2=6$ 以上の数で割り切れないのは明らかなので確かめる必要はありません)。
受け取った引数が素数であるときに True を、それ以外の数 (負数、1、合成数) を受け取ったときに False を返す素数判定関数を実装してみましょう (後述する SymPy パッケージに isprime() という関数があるので、関数名の重複を避けるためにアンダースコアを入れて is_prime() としておきます)。
# PYTHON_PRIME_NUMBER
# In[1]
# 素数判定関数(1)
def is_prime(n):
if n < 2:
return False
for k in range(2, int(n/2)+1):
if n % k == 0:
return False
return True$19$ が素数であることを確認してみましょう。
# In[2]
is_prime(19)True
合成数 $24$ に対して False を返すことを確かめておきます。
# In[3]
is_prime(24)False
より大きな数を渡してみましょう。$997$ は素数でしょうか？
# In[4]
is_prime(997)True
素数ですね。実は素数表から拾ってきた数を入れて、プログラムは正常に動くことを確認したのです。素数判定プログラムは計算効率が大切です。Jupyter notebook のマジックコマンド timeit を使って、実行時間を計ってみましょう。
%timeit is_prime(997)48.1 µs ± 3.77 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
10000 ループを 7 回走らせて、1 ループ当たりの平均実行時間が 48.1 マイクロ秒となっています。コードを少し改良するだけで、劇的に計算時間を短縮できます。初等整数論の定理によると「合成数 $N$ は $1$ より大きく $\sqrt{N}$ 以下の約数を必ずもつ」ことが知られています。
$2$ から順に $\sqrt{N}$ まで割れば、素数か合成数かを判定できるということです。これを踏まえてコードを書き直してみます。
# In[5]
# 素数判定関数(2)
def is_prime(n):
if n < 2:
return False
for k in range(2, int(n**0.5)+1):
if n % k == 0:
return False
return True先程と同じく、$997$ の素数判定の実行時間を計ってみます。
%timeit is_prime(997)3.55 µs ± 9.75 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
平均で 3.55 マイクロ秒でした。大幅に計算効率が改善されていることがわかります。is_prime() 関数を使って、指定した数以下の素数をすべて出力する関数を作成してみましょう。
# In[6]
# 素数リスト生成関数
def prime_list(m):
x = []
for k in range(m):
if is_prime(k):
x.append(k)
return x, len(x)戻り値は素数のリストと要素の数です。試しに、$50$ 以下の素数を並べてみましょう。
# In[7]
# 50以下の素数のリストを生成
prime_list(50)[0][2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
SymPy には自然数を素数であるかどうかを判定したり、素数を生成する関数がたくさん揃っています。以下に SymPy で利用できる素数関連の関数をまとめておきます。
sympy.primerange()
sympy.primerange(a, b) は a 以上 b 未満の素数を生成するジェネレータを返します。
# SYMPY_PRIME_NUMBER
# In[1]
import sympy
# 2～29の範囲内の素数を返すジェネレータ
p = sympy.primerange(2, 30)
print(p)<generator object primerange at 0x059A4510>
ジェネレータをリストに渡すと中身を列挙します。
# In[2]
# 2～29の範囲内の素数を列挙
p_list = list(p)
print(p_list)[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
sympy.prime()
sympy.prime(n) は (最小素数 2 から数えて) n 番目の素数を返します。
# In[3]
# 5番目の素数
p = sympy.prime(5)
print(p)11
sympy.nextprime()
sympy.nextprime(a, b) は a より大きな b 番目の素数を返します。
# In[4]
# 6より大きな4番目の素数
p = sympy.nextprime(6, 5)
print(p)19
sympy.prevprime()
sympy.prevprime(n) は n 以下の最大素数を返します。
# In[5]
# 30以下の最大素数
p = sympy.prevprime(30)
print(p)29
sympy.primepi()
sympy.primepi(n) は n 以下の素数の個数を返します。
# In[6]
# 30以下の素数の個数
p = sympy.primepi(30)
print(p)10
sympy.isprime()
sympy.isprime(n) は n が素数であるときは True, 素数でないとき (1 あるいは合成数であるとき) には False を返します。
# In[7]
# 29が素数であるかを判定
p = sympy.isprime(29)
print(p)True
sympy.randprime()
sympy.randprime(a, b) は a 以上 b 未満のランダムな素数を返します。
# In[8]
# 100以上200未満のランダムな素数
p = sympy.randprime(100, 200)
print(p)127
sympy.primorial()
sympy.primorial(n) は n 番目までの素数の積 (素数階乗) を返します。
# In[9]
# 5番目までの素数の積(素数階乗)
# 2*3*5*7*11
q = sympy.primorial(5)
print(q)2310
合成数と素因数分解
$1$ と素数を除く数を 合成数 (composite number) とよびます。合成数は素数の積で表せます。たとえば、$126$ は合成数であり、
\[126=2\times 3\times 3\times 7=2\times 3^2\times 7\]
のように書き表せます。合成数を構成する個々の素数を 素因数 といい、合成数を素因数の積で表すことを 素因数分解 とよびます。
素因数分解は手計算であろうと数値計算であろうと、（中学校で習ったような）「与えられた合成数を、小さい順に並べた素数で順に試し割りする」という極めて単純なアルゴリズムで実行するしかありません。もちろん、このような単純作業はコンピュータは得意としますが、与えられた合成数があまりに大きいと、それなりに時間がかかってしまうので、できるだけ効率的なプログラムを組む必要があります。
SymPy に素因数分解する関数が用意されていますが、良い練習素材なので、まずは外部ライブラリなしで素因数分解関数を実装してみましょう。素数判定プログラムに比べると、素因数分解プログラムはやや複雑です。まとめて実装しようとするとややこしいので、まず最初に「自然数 n, a およびリスト ls を受け取って、a が素数であれば n を a で試しに割ってみて、もし割り切れれば a を素因数リストに加えて再度 a で割る ... 」という処理を施す関数を用意します。つまり、素数候補を 1 個の a に限定しておいて、後で別の関数の中でこの a を動かそうということです。a の素数判定には記事前半で掲載した is_prime() 関数を使います。
# PYTHON_FACTORIZATION
# In[1]
# 素数判定関数
def is_prime(n):
if n < 2:
return False
for k in range(2, int(n**0.5)+1):
if n % k == 0:
return False
return True
# 素因数探索関数
def prime_factor(n, a, ls):
n = int(n)
a = int(a)
if is_prime(a):
while True:
# nがaで割り切れるならば
if n % a == 0:
# リストの末尾にaを追加する
ls.append(a)
# nをaで割って更新
n //= a
else:
break
return ls, n合成数 440 の中に素因数 2 が何個あるのか確かめてみましょう。
# In[2]
# 合成数440を2で繰り返し割り、
# 割り切れた回数だけリストに2を追加する
prime_factor(440, 2, [])([2, 2, 2], 55)
タプルの中のリストに素因数が記録されています。440 は 2 で 3 回割り切れること、つまり素因数 2 が 3 個含まれていることがわかりました。55 という数字は、440 を 2 で 3 回割った後の被除数です。これ以上割り切れないことを確認するために出力しています。それでは、素因数分解関数を実装してみましょう。
# In[3]
# 素因数分解関数
def factorization(n):
x = []
for k in range(2, n//2+1):
prime_factor(n, k, x)
return xfactorization() は、prime_factor() の素因数候補 a を 2 から n/2 まで動かして、それぞれの a について素数判定を行ない、素数である場合は、n を a で割り切れるごとにリストの末尾に a を付け加えます。567144 を素因数分解してみましょう。
# In[4]
# 567144を素因数分解
factorization(567144)[2, 2, 2, 3, 3, 7877]
実行結果により、$567144=2^3\times 3^2\times 7877$ の形に素因数分解できることがわかりました。
SymPy には合成数を返したり、素因数分解を実行してくれる関数があります。
sympy.composite()
sympy.composite(n) は n 番目の合成数を返します。
# SYMPY_FACTORIZATION
# In[1]
import sympy
# 5番目の合成数
c = sympy.composite(5)
print(c)10
sympy.compositepi()
sympy.compositepi(n) は n 以下の合成数の個数を返します。
# In[2]
# 12以下の合成数の個数
n = sympy.compositepi(12)
print(n)6
sympy.factorint()
sympy.factorint(n) は n を素因数分解して、結果をディクショナリで返します。
# In[3]
# 23100を素因数分解
f = sympy.factorint(23100)
print(f){2: 2, 3: 1, 5: 2, 7: 1, 11: 1}
key は素因数、value は素因数の数です。上のコードの実行結果は、$23100=2^2\times 3^1\times 5^2\times 7^1\times 11^1$ の形に素因数分解されたことを意味しています。