「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。今月(2022年10月)からは多忙につき、日々の投稿数を減らします。

MENU

Pythonで学ぶアルゴリズム(17/18)

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

を参照していく。

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

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

9.3 素数判定アルゴリズム

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

(1) 自然数kを与える。
(2) i=1とする。
(3) 2\leq x\leq\sqrt{n}を満たすような自然数xをランダムに取り、nxで割った余りが0ならば(5)に、そうでなければ(4)に移る。
(4) i\gt kならば(6)に移り、そうでなければi=i+1として(3)に戻る。
(5) n合成数であると返す。
(6) n素数であると返す。
################
### 素数判定 ###
################

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 合成数候補の探索

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

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

### 合成数かを判定するための候補がいくつあるか
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素数判定法

 \mathrm{Rabin}-\mathrm{Miller}素数判定法は如何なる合成数であっても\displaystyle{\frac{3}{4}}以上の確率で合成数かを判定するための良い候補を探すことができる方法を発見した。
 nを基数として、奇数q\leq nを用いてn-1=2^k qと表す。以下の条件

   (a) a^q\not\equiv1\ (\mathrm{mod}\ n)
   (b) {}^{\forall}i\in\{0,1,2,\cdots,k-1\}\left(a^{2^{i}q}\not\equiv-1\ (\mathrm{mod}\ n) \right)

が両方とも成立するとき、n合成数である。この判定法はここで見つけたaがどのような奇数の合成数についても75\%以上確保できることを示した点である。そのため10回連続で候補でない数を取る確率は9.5\times10^{-7}である。

###################################
### 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) # 誤り率
プライバシーポリシー お問い合わせ