素数判定法の結果をエクセルに出力したい

Question

実現したいこと

実装した3つの素数判定法の結果と計測時間をエクセルに出力したいです。フェルマーテスト、ミラーラビンテスト、aks判定法の3つを実際に実装しており、その判定結果と計測時間のプログラムも実装済みです。
理想としては、エクセルの1行目を左から、判定する数、フェルマー結果、フェルマー時間、ミラーラビンテスト結果...のようにしたいです。

発生している問題・分からないこと

現在はある一つの数について判定をするようになっていますが、そこは気にせずに教えていただきたいです。

該当のソースコード

#フェルマーテスト
import math
import random



def is_prime(n):
    if n == 1: 
        return False
    if n == 2:
        print("%d is prime." % n)
        return True

    passed_all = True

    for i in range(20):  # 20は試行回数
        a = random.randint(2, n - 1)

        if math.gcd(n, a) != 1: #互いに素ではないから合成数
            continue

        if pow(a, n - 1, n) != 1:
            passed_all = False
            #print("%d is not passed" % a)
            return False
        else:
            #print("%d is passed" % a)
            continue

    if passed_all:
        #print("%d is probably prime." % n)
        return True

    return True

#ミラーラビンテスト
import random
import math

def MillerRabinTest(n):  
    if n == 1:
        return False
    
    if n == 2:
        #print("%d is prime." % n)
        return True
    
    if n % 2 == 0:
        #print("%d is not prime" % n)
        return False
    
    passed_all = True

    s, d = CountSd(n)
    
    for i in range(20): #試行回数
        a = random.randint(2, n - 1)

        if pow(a, d, n) == 1:
            #print("%d is passed." % a)
            continue
        
        u = 0

        while u < s:
            if pow(a, d * (2 ** u), n) == n - 1: #-1=n-1
                #print("%d is passed" % a)
                break
            u += 1
        
        else:
            passed_all = False
            #print("%d is not passed" % a)
            return False
    
    if passed_all:
        #print("%d is probably prime." % n)
        return True
    
    return True
        
        
    

def CountSd(n):
    if n % 2 == 0:
        return False
    if n % 2 == 1:
        d = n - 1
        s = 0
        while d % 2 == 0:
            d //= 2
            s += 1
    return s, d
    

#aks
import math
import numpy as np

def is_prime1(n):
    if n == 1: return False

    # Step 1
    if is_perfect_power(n): return False

    # Step 2
    r = enough_order_modulo(n)

    # Step 3
    for a in range(2, min(r+1, n)):
        if n % a == 0:
            return False

    # Step 4
    if n <= r: return True

    # Step 5
    for a in range(1, int(math.sqrt(tortient(r)) * math.log(n)) + 1):
        if not is_congruent(a, n, r):
            return False

    # Step 6
    return True

# nが累乗数かどうかを判定
def is_perfect_power(n):
    size = int(math.log2(n) + 1)
    for k in range(2, size):
        pn = int(math.pow(n, 1 / k))
        if pn ** k == n or (pn+1) ** k == n: return True
    return False

# 十分大きな位数を持つ最小の数
def enough_order_modulo(n):
    a = int(math.log(n) ** 2)
    for r in range(1, n):
        order = 0
        prod = 1
        for e in range(1, r):
            prod = prod * n % r
            if prod == 1:
                order = e
                break
        if order > a:
            return r
    return n

# オイラーのトーシェント関数
def tortient(r):
    ps = primes(r)
    res = r
    for p in ps:
        res  = res * (p - 1) / p
    return res

# rの素因数のリスト
def primes(r):
    n = r
    res = set()
    for p in range(2, int(math.sqrt(r)) + 1):
        while n % p == 0:
            res.add(p)
            n /= p
    return res

# (x + a)^n と x^n + a が mod x^r - 1, n で合同かどうかを判定
def is_congruent(a, n, r):
    p = PolynomialModulo()
    ls1 = p.pow([a, 1], n, r)

    i = n % r
    ls2 = [0] * (i + 1)
    ls2[0] = a % n
    ls2[i] = 1

    return ls1 == ls2

class PolynomialModulo(object):

    def pow(self, ls, n, r):
        self.ls = ls
        self.n = n
        self.r = r

        return self.__pow(self.n)

    def __pow(self, m):
        if m == 1: return self.ls
        if m % 2 == 0:
            pls = self.__pow(m // 2)
            return self.__product(pls, pls)
        else:
            return self.__product(self.__pow(m - 1), self.__pow(1))

    def __product(self, ls1, ls2):
        res = [0] * min(len(ls1) + len(ls2) - 1, self.r)
        for i in range(len(ls1)):
            for j in range(len(ls2)):
                res[(i + j) % self.r] += ls1[i] * ls2[j]
        l = len(res)
        for k in reversed(range(l)):
            res[k] %= self.n
            if k == len(res) - 1 and res[k] == 0: res.pop(k)
        return res

import fermat_test
import MillerRabinTest
import aks_prime_test
import time


start = time.perf_counter()
print(fermat_test.is_prime(561))
end = time.perf_counter()
print(end-start)
#処理時間出力

start = time.perf_counter()
print(MillerRabinTest.MillerRabinTest(561))
end = time.perf_counter()
print(end-start)
#処理時間出力

start = time.perf_counter()
print(aks_prime_test.is_prime1(561))
end = time.perf_counter()
print(end-start)
#処理時間出力

試したこと・調べたこと

• teratailやGoogle等で検索した
• ソースコードを自分なりに変更した
• 知人に聞いた
• その他

上記の詳細・結果

データをエクセルに出力する方法を調べたが、データを配列にする必要があるとおもい、理解できず断念した。

補足

特になし

Comments

[meg_]

> データをエクセルに出力する方法を調べたが、データを配列にする必要があるとおもい、理解できず断念した。
どういうことでしょうか？実際に調べたサイトや書籍等の情報を追記いただけないでしょうか？

[jimbe]

結果を csv でファイル出力すれば EXCEL で見れると思いますが、それではダメなのでしょうか。

Answer 1

以下のコードでは openpyxl - A Python library to read/write Excel 2010 xlsx/xlsm files を利用しています。

import fermat_test
import MillerRabinTest
import aks_prime_test
import time
from openpyxl import Workbook

header = [
    '判定する数',
    'フェルマー結果', 'フェルマー時間',
    'ミラーラビンテスト結果', 'ミラーラビンテスト時間',
    'aks判定法結果', 'aks判定法時間'
]
nums = [561, 49727, 3]
data = [header]
for n in nums:
    row = [n]
    for test_func in (fermat_test.is_prime, MillerRabinTest.MillerRabinTest, aks_prime_test.is_prime1):
        start = time.perf_counter()
        result = test_func(n)
        end = time.perf_counter()
        row += [result, end - start]
    data.append(row)

# Save to Excel spreadsheet
wb = Workbook()
ws = wb.active
for row in data:
    ws.append(row)

ws.freeze_panes = 'A2'
wb.save('primality_test_result.xlsx')