Pythonで素数を扱うクラスを作る

はじめに

paizaラーニング　レベルアップ問題集「素数メニュー」の「素数大学」を解いたので、解法をクラスにまとめました。
問題とは直接関係のないメソッドも入れてあります。

プログラミング練習問題「グロタンディーク素数」

素因数分解は次の「最小公倍数」の範囲ですが、取り組む前になんとなく書けてしまったのでこのクラスに入れました。

クラス定義

変数

• NMAX: 問題集で扱う数の最大数（600万）

• isprime[]: 素数判定リスト

• primes[]: 素数のみのリスト

メソッド

• is_prime(n): 素数判定

• num_primes(): オブジェクトで扱える素数の数

• nth_prime(n): n番目の素数

• largest_prime(): オブジェクトで扱える最大の素数

• prime_factorize(n): nを素因数分解する

• is_probably_prime(n): フェルマーテストを用いて確率的素数判定をする

import random

class MyPrimeNumber:
    def __init__(self):
        self.NMAX = 6 * 10**6

        # 素数判定リストを作る
        self.isprime = [True] * (self.NMAX + 1)
        self.set_isprime()

        # 素数リストを作る
        self.primes = []
        self.genprimes()

    def set_isprime(self):
        # 素数判定リスト作成
        self.isprime[0] = False
        self.isprime[1] = False
        for i in range(2, self.NMAX + 1):
            if self.isprime[i]:
                k = 2
                while i * k <= self.NMAX:
                    self.isprime[i * k] = False
                    k += 1
                i += 1

    def is_prime(self, n: int):
        # 素数判定
        return self.isprime[n]

    def genprimes(self):
        # 素数リスト作成
        for num in range(self.NMAX + 1):
            if self.is_prime(num):
                self.primes.append(num)

    def num_primes(self):
        # このクラスで作られる素数の数を返す
        return len(self.primes)

    def nth_prime(self, n: int):
        # n番目の素数を返す
        return self.primes[n - 1]

    def largest_prime(self):
        # このクラスで扱える最大の素数を返す
        return self.primes[-1]

    def prime_factorize(self, n: int):
        # nを素因数分解しリストで返す
        if n <= 1:
            return None
        elif n <= 3:
            return [n]
        for number in range(4, n + 1):
            factors = []
            i = 0
            while self.primes[i] <= number:
                if number % self.primes[i] == 0:
                    factors.append(self.primes[i])
                    number //= self.primes[i]
                else:
                    i += 1
        return factors

    def power(self, a, n, p):
        # (a^n) % pを効率良く求める
        res = 1
        a = a % p

        while n > 0:
            if n % 2:
                res = (res * a) % p
                n -= 1
            else:
                a = a**2 % p
                n //= 2

        return res % p

    def is_probably_prime(self, n: int):
        # フェルマーテスト
        # nが確率的素数であるか判定
        if n < 2:
            return False
        if n == 2:
            return True

        # 2以上n-1以下の乱数aでnが割り切れれば合成数
        a = random.randint(2, n)
        if n % a == 0:
            return False

        # a^(n-1) % n が1であれば「たぶん素数」と判定
        fermat = self.power(a, n-1, n)
        if fermat == 1:
            return True
        else:
            return False

使用例

オブジェクトを作成します。

p = MyPrimeNumber()

素数判定

57はグロタンディーク素数と呼ばれる数です。

print(p.is_prime(19))
print(p.is_prime(57))    # グロタンディーク素数（合成数）

True
False

素数リストのサイズ、n番目の素数、オブジェクトで最大の素数

print(p.num_primes())    # 素数の数
print(p.nth_prime(50))   # 50番目の素数
print(p.largest_prime()) # 最大の素数

412849
229
5999993

素因数分解

print(p.prime_factorize(72))
print(p.prime_factorize(1))    # 素因数分解できない
print(p.prime_factorize(3))    # 素数には単一要素のリストを返す
print(p.prime_factorize(5963))

[2, 2, 2, 3, 3]
None
[3]
[67, 89]

フェルマーテスト（確率的素数判定）

問題集では底を2に固定していましたが、底に乱数を採用する手法に変えました。

91をフェルマーテストにかけてみます。

k = 91

n = 100
counter = 0
for _ in range(n):
    if p.is_probably_prime(k):
        counter += 1
print(counter/n)    # 100回中何回素数と判定されるかで確率を出す

print(p.prime_factorize(k))

素数の確率が半々程度です。実際に7×13の合成数です。

0.41
[7, 13]

2821はフェルマーテストにおいて高い確率で素数と判定されますが、合成数です。

k = 2821

n = 100
counter = 0
for _ in range(n):
    if p.is_probably_prime(k):
        counter += 1
print(counter/n)

print(p.prime_factorize(k))

0.84
[7, 13, 31]

おわりに

ローカルで組むなら、もっと大きい数までの素数を求めて外部テキストファイルに書き出し、素数判定はファイル読み込みを使うのが良いだろうと思います。

フェルマーテストのほかに「ミラー-ラビン素数判定法」などの素数判定法があるので作ってみようかと思いましたが、問題集の範囲にはなさそうなので今後の楽しみにとっておきます。