結果

問題 No.981 一般冪乗根
ユーザー hari64
提出日時 2021-07-22 15:24:31
言語 PyPy3
(7.3.15)
結果
AC  
実行時間 3,656 ms / 6,000 ms
コード長 6,783 bytes
コンパイル時間 359 ms
コンパイル使用メモリ 82,576 KB
実行使用メモリ 215,412 KB
最終ジャッジ日時 2024-07-17 14:53:47
合計ジャッジ時間 164,466 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 30 TLE * 2 MLE * 12
権限があれば一括ダウンロードができます

ソースコード

diff # 

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a


def isPrimeMR(n):
    d = n - 1
    d = d // (d & -d)
    L = [2]
    for a in L:
        t = d
        y = pow(a, t, n)
        if y == 1:
            continue
        while y != n - 1:
            y = (y * y) % n
            if y == 1 or t == n - 1:
                return 0
            t <<= 1
    return 1


def findFactorRho(n):
    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 isPrimeMR(g):
                return g
            elif isPrimeMR(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 isPrimeMR(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 prime_factors(n, return_e: bool = True):
#     """
#     nの素因数分解を指数付きで返す 計算量O(√N) O(n^(1/4))もある
#     https://qiita.com/drken/items/a14e9af0ca2d857dad23#4-%E7%B4%A0%E5%9B%A0%E6%95%B0%E5%88%86%E8%A7%A3
#     """
#     assert 1 <= n
#     pfs = []
#     a = 2
#     while a**2 <= n:
#         if n % a == 0:
#             ex = 0
#             while n % a == 0:
#                 ex += 1
#                 n //= a
#             pfs.append((a, ex) if return_e else a)
#         a += 1
#     if n != 1:
#         pfs.append((n, 1) if return_e else n)
#     return pfs


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