結果

問題 No.3030 ミラー・ラビン素数判定法のテスト
ユーザー JashinchanJashinchan
提出日時 2022-09-02 12:06:55
言語 Python3
(3.12.2 + numpy 1.26.4 + scipy 1.12.0)
結果
RE  
実行時間 -
コード長 11,699 bytes
コンパイル時間 81 ms
コンパイル使用メモリ 14,464 KB
実行使用メモリ 12,800 KB
最終ジャッジ日時 2024-04-27 14:47:43
合計ジャッジ時間 11,785 ms
ジャッジサーバーID
（参考情報） 		judge3 / judge5
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 32 ms
12,672 KB
testcase_01 AC 34 ms
12,672 KB
testcase_02 AC 33 ms
12,672 KB
testcase_03 AC 32 ms
12,672 KB
testcase_04 RE -
testcase_05 RE -
testcase_06 AC 1,841 ms
12,800 KB
testcase_07 AC 1,846 ms
12,672 KB
testcase_08 AC 1,814 ms
12,672 KB
testcase_09 RE -
権限があれば一括ダウンロードができます

ソースコード

diff # 

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)))
0