Pythonで学ぶアルゴリズム（17/18）

　アルゴリズムを学ぶのにPythonの学習も兼ねて

Ｐｙｔｈｏｎで学ぶアルゴリズムとデータ構造 (データサイエンス入門シリーズ)

を参照していく。

前回

https://power-of-awareness.com/entry/2022/07/18/120000power-of-awareness.com

9.　乱択アルゴリズムと数論

　実行の一部にランダムな挙動を含むアルゴリズムを乱択アルゴリズムという。乱択アルゴリズムを扱うと効率的な解法をもたらすような問題が存在する。

9.3　素数判定アルゴリズム

　自然数が素数かを判定するモンテカルロ・アルゴリズムを考えてみる。

(1)	自然数を与える。
(2)	とする。
(3)	を満たすような自然数をランダムに取り、をで割った余りがならば(5)に、そうでなければ(4)に移る。
(4)	ならば(6)に移り、そうでなければとして(3)に戻る。
(5)	は合成数であると返す。
(6)	は素数であると返す。

################
### 素数判定 ###
################

import random

def random_div(n, repeat = 10):
    """与えられた数nが素数かどうかを判定する。
    ランダムに割り算をして約数を見つける
    """
    
    if n % 2 ==0:
        return 'composite'
    d_max = int(pow(n, 0.5))
    
    # 奇数列を作る
    odd_seq = range(3, d_max + 1, 2)
    
    for cnt in range(repeat):
        d = random.choice(odd_seq)
        
        if n % d ==0:
            return 'composite'
    return 'probably prime'

# 簡単な素数は判定できる
random_div(71)

# 4桁の自然数になるともう怪しい
for i in range(20):
    print(random_div(1105))

9.3.1　合成数候補の探索

　先ほどのアルゴリズムを考える。探索対象となる自然数を素因数分解したときの値で割ることができれば、探索対象が合同数であることは確定する。すなわち先程は割る数をランダムに取ったが、その候補に素因数が含まれる割合が多ければ、より効率よく素数か否かを判定できる。そこでの小定理を活用してみよう。
　が素数ならば、であるようなに対してであった。これの対偶を考えると、ならばは素数ではないことが分かる。このときのを新たに素数か否かを判定するための自然数とする。

########################
### 合成数候補の探索 ###
########################

### 合成数かを判定するための候補がいくつあるか
def fermat_check(n):
    cnt = 0
    
    for a in range(2, n):
        if pow(a, n-1, n)!=1:
            cnt += 1
    return cnt

### トライアル
print(fermat_check(71)) # 素数ならば0
print(fermat_check(99)) # 一般的な合成数
print(fermat_check(1105)) # 一般的な合成数

9.3.2　Rabin-Miller素数判定法

　-素数判定法は如何なる合成数であっても以上の確率で合成数かを判定するための良い候補を探すことができる方法を発見した。
　を基数として、奇数を用いてと表す。以下の条件

	(a)
	(b)

が両方とも成立するとき、は合成数である。この判定法はここで見つけたがどのような奇数の合成数についても以上確保できることを示した点である。そのため回連続で候補でない数を取る確率はである。

###################################
### Rabin-Miller primality test ###
###################################

### まずは奇数と2の冪乗に分解する
def squeeze_q(n):
    """ 奇数の引数nから1を引いた数を計算し、
    2の成分を搾り取った奇数と2の個数を返す
    """

    k = 0
    x = n-1
    while x % 2 !=1:
        k += 1
        x //= 2
    return (n-1)//pow(2, k), k

def rabin_miller(n, repeat = 10):
    if n<2:
        return 'give me more than 1.'
    if n==2:
        return 'prime'
    if n%2==0:
        return 'composite'
    q, k = squeeze_q(n)
    cnt = 0
    while cnt < repeat:
        a = random.randint(2, n-1)
        
        # 1つ目の条件
        cond_1 = pow(a, q, n) != 1
        temp = []
        for i in range(k):
            y = pow(2, i) * q
            c = pow(a, y, n) != n-1
            temp.append(c)
        
        # 2つ目の条件
        cond_2 = all(temp)
        
        if cond_1 and cond_2:
            return 'composite'
        cnt += 1
    return 'probably prime'

### テスト
rabin_miller(71)

num_trial = 1000
S = 0
for i in range(num_trial):
    if rabin_miller(1105)=='probably prime':
        S += 1
print(S/num_trial) # 誤り率