結果
| 問題 | No.3030 ミラー・ラビン素数判定法のテスト |
| ユーザー |
 Jashinchan
|
| 提出日時 | 2022-09-02 11:09:10 |
| 言語 | Python3 (3.12.2 + numpy 1.26.4 + scipy 1.12.0) |
| 結果 |
RE
|
| 実行時間 | - |
| コード長 | 10,298 bytes |
| コンパイル時間 | 130 ms |
| コンパイル使用メモリ | 14,080 KB |
| 実行使用メモリ | 12,288 KB |
| 最終ジャッジ日時 | 2024-04-27 14:07:51 |
| 合計ジャッジ時間 | 5,197 ms |
|
ジャッジサーバーID （参考情報） |
judge5 / judge4 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 35 ms
12,160 KB |
| testcase_01 | AC | 37 ms
12,160 KB |
| testcase_02 | AC | 34 ms
12,160 KB |
| testcase_03 | AC | 35 ms
12,160 KB |
| testcase_04 | RE | - |
| testcase_05 | RE | - |
| testcase_06 | RE | - |
| testcase_07 | AC | 2,301 ms
12,160 KB |
| testcase_08 | RE | - |
| testcase_09 | RE | - |
ソースコード
import copy
def gcd(a, b):
while a:
a, b = b % a, a
return b
def v(q, t):
ret = 0
while (t % q == 0):
ret += 1
t //= q
return ret
def is_prime_trivial_division(n):
if n < 2:
return False
elif n == 2 or n == 3:
return True
elif n % 2 == 0:
return False
else:
i = 3
while i * i <= n:
if n % i == 0:
return False
i += 2
return True
# n limit: 18446744073709551616
t = 60
et = 6814407600
q_list = [2, 3, 5, 7, 11, 13, 31, 61]
factors_q1_list = {
2: [],
3: [(2, 1)],
5: [(2, 2)],
7: [(2, 1), (3, 1)],
11: [(2, 1), (5, 1)],
13: [(2, 2), (3, 1)],
31: [(2, 1), (3, 1), (5, 1)],
61: [(2, 2), (3, 1), (5, 1)],
}
fq_list = {
5: {1: 2, 2: 1, 3: 3},
7: {1: 5, 2: 3, 3: 2, 4: 4, 5: 1},
11: {1: 5, 2: 3, 3: 2, 4: 7, 5: 1, 6: 8, 7: 4, 8: 6, 9: 9},
13: {1: 6, 2: 10, 3: 5, 4: 7, 5: 3, 6: 1, 7: 4, 8: 9, 9: 8, 10: 2, 11: 11},
31: {1: 9, 2: 27, 3: 20, 4: 11, 5: 25, 6: 6, 7: 21, 8: 19, 9: 1, 10: 28, 11: 4, 12: 13, 13: 12, 14: 17, 15: 24, 16: 18, 17: 14, 18: 16, 19: 8, 20: 3, 21: 7, 22: 26, 23: 29, 24: 15, 25: 5, 26: 22, 27: 2, 28: 10, 29: 23},
61: {1: 30, 2: 36, 3: 19, 4: 58, 5: 29, 6: 31, 7: 52, 8: 45, 9: 27, 10: 50, 11: 18, 12: 33, 13: 17, 14: 41, 15: 53, 16: 25, 17: 13, 18: 11, 19: 3, 20: 28, 21: 35, 22: 32, 23: 42, 24: 56, 25: 16, 26: 43, 27: 9, 28: 20, 29: 5, 30: 1, 31: 6, 32: 22, 33: 12, 34: 47, 35: 21, 36: 2, 37: 49, 38: 40, 39: 44, 40: 38, 41: 14, 42: 23, 43: 26, 44: 39, 45: 8, 46: 57, 47: 34, 48: 51, 49: 37, 50: 10, 51: 48, 52: 7, 53: 15, 54: 55, 55: 54, 56: 24, 57: 46, 58: 4, 59: 59},
}
class JacobiSum(object):
def __init__(self, p, k, q):
self.p = p
self.k = k
self.q = q
self.m = (p - 1) * p**(k - 1)
self.pk = p**k
self.coef = [0] * self.m
def one(self):
self.coef[0] = 1
for i in range(1, self.m):
self.coef[i] = 0
return self
def mul(self, jac):
m = self.m
pk = self.pk
j_ret = JacobiSum(self.p, self.k, self.q)
for i in range(m):
for j in range(m):
if (i + j) % pk < m:
j_ret.coef[(i + j) % pk] += self.coef[i] * jac.coef[j]
else:
r = (i + j) % pk - self.p**(self.k - 1)
while r >= 0:
j_ret.coef[r] -= self.coef[i] * jac.coef[j]
r -= self.p ** (self.k - 1)
return j_ret
def __mul__(self, right):
if type(right) is int:
j_ret = JacobiSum(self.p, self.k, self.q)
for i in range(self.m):
j_ret.coef[i] = self.coef[i] * right
return j_ret
else:
return self.mul(right)
def modpow(self, x, n):
j_ret = JacobiSum(self.p, self.k, self.q)
j_ret.coef[0] = 1
j_a = copy.deepcopy(self)
while x > 0:
if x % 2 == 1:
j_ret = (j_ret * j_a).mod(n)
j_a = j_a * j_a
j_a.mod(n)
x //= 2
return j_ret
def mod(self, n):
for i in range(self.m):
self.coef[i] %= n
return self
def sigma(self, x):
m = self.m
pk = self.pk
j_ret = JacobiSum(self.p, self.k, self.q)
for i in range(m):
if (i * x) % pk < m:
j_ret.coef[(i * x) % pk] += self.coef[i]
else:
r = (i * x) % pk - self.p ** (self.k - 1)
while r >= 0:
j_ret.coef[r] -= self.coef[i]
r -= self.p ** (self.k - 1)
return j_ret
def sigma_inv(self, x):
m = self.m
pk = self.pk
j_ret = JacobiSum(self.p, self.k, self.q)
for i in range(pk):
if i < m:
if (i * x) % pk < m:
j_ret.coef[i] += self.coef[(i * x) % pk]
else:
r = i - self.p ** (self.k - 1)
while r >= 0:
if (i * x) % pk < m:
j_ret.coef[r] -= self.coef[(i * x) % pk]
r -= self.p ** (self.k - 1)
return j_ret
def is_root_of_unity(self, N):
m = self.m
p = self.p
k = self.k
one = 0
for i in range(m):
if self.coef[i] == 1:
one += 1
h = i
elif self.coef[i] == 0:
continue
elif (self.coef[i] - (-1)) % N != 0:
return False, None
if one == 1:
return True, h
for i in range(m):
if self.coef[i] != 0:
break
r = i % (p**(k - 1))
for i in range(m):
if i % (p**(k - 1)) == r:
if (self.coef[i] - (-1)) % N != 0:
return False, None
else:
if self.coef[i] != 0:
return False, None
return True, (p - 1) * p**(k - 1) + r
def calc_J_ab(p, k, q, a, b):
j_ret = JacobiSum(p, k, q)
f = fq_list[q]
for x in range(1, q - 1):
pk = p**k
if (a * x + b * f[x]) % pk < j_ret.m:
j_ret.coef[(a * x + b * f[x]) % pk] += 1
else:
r = (a * x + b * f[x]) % pk - p**(k - 1)
while r >= 0:
j_ret.coef[r] -= 1
r -= p**(k - 1)
return j_ret
def calc_J(p, k, q):
return calc_J_ab(p, k, q, 1, 1)
def calc_J3(p, k, q):
j2q = calc_J(p, k, q)
j21 = calc_J_ab(p, k, q, 2, 1)
j_ret = j2q * j21
return j_ret
def calc_J2(p, k, q):
j31 = calc_J_ab(2, 3, q, 3, 1)
j_conv = JacobiSum(p, k, q)
for i in range(j31.m):
j_conv.coef[i * (p**k) // 8] = j31.coef[i]
j_ret = j_conv * j_conv
return j_ret
def step4a(p, k, q, n):
J = calc_J(p, k, q)
s1 = JacobiSum(p, k, q).one()
for x in range(p**k):
if x % p == 0:
continue
t = J.sigma_inv(x)
t = t.modpow(x, n)
s1 = s1 * t
s1.mod(n)
r = n % (p**k)
s2 = s1.modpow(n // (p**k), n)
J_alpha = JacobiSum(p, k, q).one()
for x in range(p**k):
if x % p == 0:
continue
t = J.sigma_inv(x)
t = t.modpow((r * x) // (p**k), n)
J_alpha = J_alpha * t
J_alpha.mod(n)
S = (s2 * J_alpha).mod(n)
exist, h = S.is_root_of_unity(n)
if not exist:
return False, None
else:
if h % p != 0:
l_p = 1
else:
l_p = 0
return True, l_p
def step4b(p, k, q, n):
J = calc_J3(p, k, q)
s1 = JacobiSum(p, k, q).one()
for x in range(p**k):
if x % 8 not in [1, 3]:
continue
t = J.sigma_inv(x)
t = t.modpow(x, n)
s1 = s1 * t
s1.mod(n)
r = n % (p**k)
s2 = s1.modpow(n // (p**k), n)
J_alpha = JacobiSum(p, k, q).one()
for x in range(p**k):
if x % 8 not in [1, 3]:
continue
t = J.sigma_inv(x)
t = t.modpow((r * x) // (p**k), n)
J_alpha = J_alpha * t
J_alpha.mod(n)
if n % 8 in [1, 3]:
S = (s2 * J_alpha).mod(n)
else:
J2_delta = calc_J2(p, k, q)
S = (s2 * J_alpha * J2_delta).mod(n)
exist, h = S.is_root_of_unity(n)
if not exist:
return False, None
else:
if h % p != 0 and (pow(q, (n - 1) // 2, n) + 1) % n == 0:
l_p = 1
else:
l_p = 0
return True, l_p
def step4c(p, k, q, n):
J2q = calc_J(p, k, q)
s1 = (J2q * J2q * q).mod(n)
s2 = s1.modpow(n // 4, n)
if n % 4 == 1:
S = s2
elif n % 4 == 3:
S = (s2 * J2q * J2q).mod(n)
else:
exit(1)
exist, h = S.is_root_of_unity(n)
if not exist:
return False, None
else:
if h % p != 0 and (pow(q, (n - 1) // 2, n) + 1) % n == 0:
l_p = 1
else:
l_p = 0
return True, l_p
def step4d(p, k, q, n):
S2q = pow(-q, (n - 1) // 2, n)
if (S2q - 1) % n != 0 and (S2q + 1) % n != 0:
return False, None
else:
if (S2q + 1) % n == 0 and (n - 1) % 4 == 0:
l_p = 1
else:
l_p = 0
return True, l_p
def step4(p, k, q, N):
if p >= 3:
result, l_p = step4a(p, k, q, N)
elif p == 2 and k >= 3:
result, l_p = step4b(p, k, q, N)
elif p == 2 and k == 2:
result, l_p = step4c(p, k, q, N)
elif p == 2 and k == 1:
result, l_p = step4d(p, k, q, N)
else:
exit(1)
return result, l_p
def APR_CL_primality_test(n):
# step 0
if n < 2:
return False
if n < 4:
return True
if n & 1 == 0:
return False
# step 1
g = gcd(408864456000, n)
if g > 1:
return False
l = {}
# step 2
factor_t = {2: 2, 3: 1, 5: 1}
for p, k in factor_t.items():
if p >= 3 and pow(n, p - 1, p * p) != 1:
l[p] = 1
else:
l[p] = 0
# step 3
for q, fac in factors_q1_list.items():
if q == 2:
continue
for p, k in fac:
# step 4
result, l_p = step4(p, k, q, n)
if not result:
return False
elif l_p == 1:
l[p] = 1
# step 5
for p, value in l.items():
if value == 0:
count = 0
i = 1
found = False
while count < 5:
q = p * i + 1
if n % q != 0 and is_prime_trivial_division(q) and (q not in q_list):
count += 1
k = v(p, q - 1)
result, l_p = step4(p, k, q, n)
if not result:
return False
elif l_p == 1:
found = True
break
i += 1
if not found:
exit(1)
# step 6
r = 1
for _ in range(59):
r = (r * n) % et
if r != 1 and r != n and n % r == 0:
return False
return True
for i in range(int(input())):
x = int(input())
print(x, int(APR_CL_primality_test(x)))