結果
| 問題 | No.3030 ミラー・ラビン素数判定法のテスト |
| ユーザー |
 Jashinchan
|
| 提出日時 | 2022-09-02 12:06:55 |
| 言語 | Python3 (3.12.2 + numpy 1.26.4 + scipy 1.12.0) |
| 結果 |
RE
|
| 実行時間 | - |
| コード長 | 11,699 bytes |
| コンパイル時間 | 328 ms |
| コンパイル使用メモリ | 14,464 KB |
| 実行使用メモリ | 12,672 KB |
| 最終ジャッジ日時 | 2024-11-15 17:35:55 |
| 合計ジャッジ時間 | 13,074 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge5 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 38 ms
12,544 KB |
| testcase_01 | AC | 41 ms
12,672 KB |
| testcase_02 | AC | 38 ms
12,544 KB |
| testcase_03 | AC | 39 ms
12,672 KB |
| testcase_04 | RE | - |
| testcase_05 | RE | - |
| testcase_06 | AC | 2,060 ms
12,672 KB |
| testcase_07 | AC | 2,032 ms
12,544 KB |
| testcase_08 | AC | 2,027 ms
12,672 KB |
| 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},
17: {1: 6, 2: 2, 3: 10, 4: 5, 5: 4, 6: 1, 7: 11, 8: 14, 9: 12, 10: 3, 11: 7, 12: 9, 13: 15, 14: 8, 15: 13},
29: {1: 14, 2: 19, 3: 26, 4: 13, 5: 15, 6: 8, 7: 11, 8: 6, 9: 25, 10: 17, 11: 7, 12: 20, 13: 4, 14: 1, 15: 5, 16: 22, 17: 10, 18: 21, 19: 2, 20: 12, 21: 18, 22: 16, 23: 24, 24: 23, 25: 9, 26: 3, 27: 27},
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},
37: {1: 18, 2: 8, 3: 14, 4: 31, 5: 27, 6: 30, 7: 22, 8: 2, 9: 32, 10: 11, 11: 10, 12: 28, 13: 15, 14: 3, 15: 13, 16: 21, 17: 25, 18: 1, 19: 26, 20: 23, 21: 16, 22: 7, 23: 20, 24: 34, 25: 17, 26: 19, 27: 5, 28: 12, 29: 33, 30: 6, 31: 4, 32: 9, 33: 29, 34: 24, 35: 35},
41: {1: 2, 2: 1, 3: 28, 4: 33, 5: 37, 6: 15, 7: 31, 8: 10, 9: 36, 10: 8, 11: 25, 12: 35, 13: 16, 14: 14, 15: 6, 16: 13, 17: 24, 18: 30, 19: 38, 20: 26, 21: 39, 22: 32, 23: 27, 24: 17, 25: 11, 26: 20, 27: 23, 28: 3, 29: 34, 30: 18, 31: 7, 32: 22, 33: 4, 34: 29, 35: 12, 36: 9, 37: 5, 38: 19, 39: 21},
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 smallest_primitive_root(q):
for r in range(2, q):
s = set({})
m = 1
for i in range(1, q):
m = (m * r) % q
s.add(m)
if len(s) == q - 1:
return r
return None
def calc_f(q):
g = smallest_primitive_root(q)
m = {}
for x in range(1, q - 1):
m[pow(g, x, q)] = x
f = {}
for x in range(1, q - 1):
f[x] = m[(1 - pow(g, x, q)) % q]
return f
def calc_J_ab(p, k, q, a, b):
j_ret = JacobiSum(p, k, q)
if q in q_list:
f = fq_list[q]
else:
f = calc_f(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.items():
# 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)))