prime_test.py

#!/usr/bin/env python3

import timeit
import random

# ここでは確率的素数判定法を実装する
# 確率的素数判定法は合成数かそうではなさそうか判断する。複数回行って、信用性を確保するみたい

# まずは、原始的な考えで
def way1(num):
    """numが素数かどうかチェックする。
    順に割っていくので桁が大きくなると、とても遅い
    """
    st = timeit.default_timer()
    # 2からnum-1まで順に割っていき、もし、割りきれたら、素数でない
    isPrime = True
    aFactor = 0
    for i in range(2, num-1):
        if num % i == 0:
            isPrime = False
            aFactor = i
            break
    end = timeit.default_timer()
    # if isPrime:
    #     print(num, 'is prime.', end - st, 'ms')
    # else:
    #     print(num, 'is composite. A factor is', aFactor, ',', end - st, 'ms')
    return isPrime
    # print('Proccess time is', end - st, 'ms')

def way2(n, k=20):
    """n > 2が素数か判定する
    ミラー–ラビン素数判定法（Miller–Rabin primality test）
    """
    li = []
    if n % 2 == 0:
        return False
    s = 0
    n1 = n - 1
    d = n1
    # Ex. n=9, n-1=8, 8=2**3*1
    # Ex. n=15, n-1=14, 14=2**1*7 ?
    # Ex. n=27, n-1=26, 26=2**1*13 ?
    while d % 2 == 0:
        s += 1
        d //= 2
    # n-1=2**s*d
    # print(n1,'= ( 2 ^', s, ') *', d)
    for i in range(k):
        a = random.randint(1, n1)
        first = 1 != pow(a, d, n)
        li.clear()
        for r in range(0, s): # 0~s-1の範囲
            result = pow(a, pow(2,r)*d, n)
            li.append(result)
        if first and li.count(n1) == 0: # n-1 = -1 mod n
            # print(n, 'is composite. Time:', end-st, 'ms, Witness:', a)
            return False # 合成数
    # print(n, 'is probably prime. Time:', end-st, 'ms')
    # 素数
    return True

    
if __name__ == '__main__':
    # for i in range(2000, 2100):
    #     way1(i)
    # print('way1')
    # way1(1001)
    # # way1(150094669656737489) # 時間がかかりすぎる15分以上？
    # way1(181243) # 素数. 0.02315494802314788 ms -> 0.023ms
    # way1(539633) # 素数. 0.08834185401792638 ms -> 0.088ms
    # print('way2')
    # way2(1001)
    # way2(150094669656737489) # ミラーラビンだと0.0008150929934345186 msで済む
    # way2(181243)
    # way2(539633)
    loopCount = 10000
    
    com = 0
    pri = 0
    com2 = 0
    pri2 = 0
    
    st = timeit.default_timer()
    for i in range(3, loopCount):
        if way1(i):
            # print(i, 'is probably prime.')
            pri += 1
        else:
            # print(i, 'is composite.')
            com += 1
    end = timeit.default_timer()
    for i in range(3, loopCount):
        if way2(i):
            pri2 += 1
        else:
            com2 += 1
    end2 = timeit.default_timer()
    print('Way1. Primes:', pri, ', Composites:', com, 'Time:', end-st, 's')
    print('Way2. Primes:', pri2, ', Composites:', com2, 'Time:', end2-end, 's')
    # 500000まで素数を探した時
    # Way1. Primes: 41537 , Composites: 458460 Time: 1196.9471503770037 s
    # Way2. Primes: 41537 , Composites: 458460 Time: 12.93861346898484 s