結果

問題 No.981 一般冪乗根
ユーザー hari64
提出日時 2021-07-22 18:09:46
言語 PyPy3
(7.3.15)
結果
AC  
実行時間 3,642 ms / 6,000 ms
コード長 6,784 bytes
コンパイル時間 1,106 ms
コンパイル使用メモリ 82,604 KB
実行使用メモリ 214,184 KB
最終ジャッジ日時 2024-10-10 07:56:23
合計ジャッジ時間 177,435 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 30 TLE * 5 MLE * 9
権限があれば一括ダウンロードができます

ソースコード

diff # 

def is_prime(n: int) -> bool: #ミラー–ラビン素数判定法
    # https://qiita.com/srtk86/items/609737d50c9ef5f5dc59
    # https://en.wikipedia.org/wiki/Strong_pseudoprime
    assert isinstance(n, int) and 0 < n
    if n == 2:
        return True
    if n == 1 or n & 1 == 0:
        return False
    d = (n - 1) >> 1
    while d & 1 == 0:
        d >>= 1  # n-1=(2**s)*d (dは奇数)
    L=[a for a in (2, 325, 9375, 28178,450775, 9780504, 1795265022) if a<n]
    if n>=2**64:  # これはwikiによる
        from random import randint
        L += [randint(1, n-1) for _ in range(100)]
    for a in L:  # nが素数ならばa^d≡1 (mod p) もしくは
        t = d  # a^((2^r)*d)≡-1 (mod p) が成立すべき
        y = pow(a, t, n)
        while t != n - 1 and y != 1 and y != n - 1:
            y = (y * y) % n
            t <<= 1
        if y != n - 1 and t & 1 == 0:
            return False
    else:
        return True


def findFactorRho(n):
    from math import gcd
    m = 1 << n.bit_length() // 8
    for c in range(1, 99):
        def f(x):
            return (x * x + c) % n
        y, r, q, g = 2, 1, 1, 1
        while g == 1:
            x = y
            for i in range(r):
                y = f(y)
            k = 0
            while k < r and g == 1:
                ys = y
                for i in range(min(m, r - k)):
                    y = f(y)
                    q = q * abs(x - y) % n
                g = gcd(q, n)
                k += m
            r <<= 1
        if g == n:
            g = 1
            while g == 1:
                ys = f(ys)
                g = gcd(abs(x - ys), n)
        if g < n:
            if is_prime(g):
                return g
            elif is_prime(n // g):
                return n // g
            return findFactorRho(g)


def prime_factors(n, return_e: bool = False):
    i = 2
    ret = {}
    rhoFlg = 0
    while i*i <= n:
        k = 0
        while n % i == 0:
            n //= i
            k += 1
        if k:
            ret[i] = k
        i += 1 + i % 2
        if i == 101 and n >= 2 ** 20:
            while n > 1:
                if is_prime(n):
                    ret[n], n = 1, 1
                else:
                    rhoFlg = 1
                    j = findFactorRho(n)
                    k = 0
                    while n % j == 0:
                        n //= j
                        k += 1
                    ret[j] = k
    if n > 1:
        ret[n] = 1
    if rhoFlg:
        ret = {x: ret[x] for x in sorted(ret)}
    return list(ret.keys())


def primitive_root(p: int) -> int:
    """素数pに対する何らかの原始根rを一つ求める (rの位数はp-1)
    https://manabitimes.jp/math/842
    https://37zigen.com/primitive-root/ 最速ではない
    """
    if p == 2:
        return 1
    from random import randint
    p_1_pfs = prime_factors(p-1, return_e=False)
    cnt = 0
    while cnt < 2*p:
        cnt += 1
        cfr = randint(2, p-1)  # candidates for r
        for pf in p_1_pfs:
            if pow(cfr, (p-1)//pf, p) == 1:  # p-1//pfが位数候補
                break
        else:  # 全ての位数候補に対して、その否定がなされた。
            return cfr
    else:
        raise AssertionError


def discrete_logarithm(x, y, mod):
    """離散対数問題 x^k≡yとなるようなkを求める O(√mod)"""
    babys = {}
    sqrt = int(mod**0.5)+1
    # baby-step
    x_baby = 1
    for i in range(sqrt):
        babys[x_baby] = i  # x^(sqrt-1)まで求めておく
        if x_baby == y:
            return i
        x_baby = x_baby * x % mod
    # giant-step
    x_giant = pow(x_baby, mod-2, mod)  # x^(-sqrt)
    for i in range(1, sqrt+1):
        y = y * x_giant % mod
        if y in babys:  # ここで√modずつ見ていく為に速い
            return babys[y] + i*sqrt
    return -1  # 原始根でない限り、必ず存在するとは言えない


def _inv_gcd(a, b) -> tuple:
    """
    返り値は(gcd(a,b),x) (ただしxa≡g (mod b) 0<=x<b//g)
    計算量 O(log(max(a,b))) ラメの定理などからも分かる通り計算量が非常に小さい
    (0<=a<=b∊N len(str(a))=dの時、gcd算出の為の計算回数は5d以下)
    逆元を求めるならpow(g,-1,b)でもいいが、versionの問題とgとbが互いに素という制約はある
    また、このコードは拡張ユークリッドの互除法の行列表現からも理解可能 再帰でも書ける
    https://github.com/atcoder/ac-library/blob/master/atcoder/internal_math.hpp
    """
    # assert 0 <= a and 1 <= b
    a %= b
    if a == 0:
        return (b, 0)
    s = b  # m0*a≡s (mod b) これらの性質をm0,m1は満たしていることに注意
    t = a  # m1*a≡t (mod b)
    m0 = 0  # s*|m1|+t*|m0|<=b
    m1 = 1  # また、a<bよりt<sが成立
    while t:
        # この三行は互除法そのもの
        u = s // t
        s -= t * u
        s, t = t, s
        # m0とm1が先ほど挙げた性質を保持し続ける様に変更する oがoldm、nがnewを示すとして
        # o_m0*a=o_s o_m1*a=o_t (冒頭にあげた性質)
        # n_s=o_t n_t=o_s-o_t*u (互除法による変更)
        # o_m1*a=n_s (o_m0-o_m1*u)*a=n_t (代入して整理)
        # これを先の性質の式と見比べるとn_m0=o_m1 n_m1=o_m0-o_m1*uを得る これが下の式の正体
        # 三つ目の性質も保たれていることは代入すればすぐわかる
        m0 -= m1 * u
        m0, m1 = m1, m0
    if m0 < 0:  # 三つ目の性質を利用するとu*|n_m0|<=b//s=b//g⇒|n_mo|<b/g (∵u>=0)
        m0 += b // s
    return (s, m0)


def eea(a: int, b: int) -> tuple:
    """
    ax+by=gcd(a,b)なる(x,y)を求める(拡張されたユークリッドの互除法)
    (extended_euclidean_algorithm)(元のユークリッドの互除法はgcdを求める手法を指す)
    https://qiita.com/drken/items/b97ff231e43bce50199a
    https://ja.wikipedia.org/wiki/ユークリッドの互除法 <-英語版が優秀すぎ
    https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
    また、以上の話題は連分数展開とも関連があるらしい
    """
    g, x = _inv_gcd(a, abs(b))
    y = (b//abs(b))*(g-a*x)//abs(b)
    assert (g-a*x) % abs(b) == 0
    return (x, y)


def main():
    import sys
    from math import gcd
    input = sys.stdin.buffer.readline
    T = int(input())
    for _ in range(T):
        p, k, a = map(int, input().split())
        g = primitive_root(p)
        y = discrete_logarithm(g, a, p)
        gkp = gcd(k, p-1)
        if y % gkp != 0:
            print(-1)
            continue
        z = (y//gkp)*eea(k//gkp, -(p-1)//gkp)[0]
        x = pow(g, z % (p-1), p)
        print(x)


if __name__ == '__main__':
    main()
0