#680969 (Python3) No.981 一般冪乗根

• 問題: No.981 一般冪乗根
• ユーザー: hari64
• 提出日時: 2021-07-22 15:11:37
• 言語: Python3 (3.12.2 + numpy 1.26.4 + scipy 1.12.0)
• 結果: RE
• 実行時間: -
• コード長: 4,836 bytes
• コンパイル時間: 305 ms
• コンパイル使用メモリ: 11,028 KB
• 実行使用メモリ: 60,816 KB
• 最終ジャッジ日時: 2023-09-24 14:02:56
• 合計ジャッジ時間: 12,622 ms
• ジャッジサーバーID （参考情報）: judge14 / judge11

入力	結果	実行時間; 実行使用メモリ
testcase_00	RE	-
testcase_01	RE	-
testcase_02	RE	-
testcase_03	RE	-
testcase_04	RE	-
testcase_05	TLE	-
testcase_06	--	-
testcase_07	--	-
testcase_08	--	-
testcase_09	--	-
testcase_10	--	-
testcase_11	--	-
testcase_12	--	-
testcase_13	--	-
testcase_14	--	-
testcase_15	--	-
testcase_16	--	-
testcase_17	--	-
testcase_18	--	-
testcase_19	--	-
testcase_20	--	-
testcase_21	--	-
testcase_22	--	-
testcase_23	--	-
testcase_24	--	-
testcase_25	--	-
testcase_26	--	-
testcase_27	--	-
testcase_28	--	-
evil_60bit1.txt	--	-
evil_60bit2.txt	--	-
evil_60bit3.txt	--	-
evil_hack	--	-
evil_hard_random	--	-
evil_hard_safeprime.txt	--	-
evil_hard_tonelli0	--	-
evil_hard_tonelli1	--	-
evil_hard_tonelli2	--	-
evil_hard_tonelli3	--	-
evil_sefeprime1.txt	--	-
evil_sefeprime2.txt	--	-
evil_sefeprime3.txt	--	-
evil_tonelli1.txt	--	-
evil_tonelli2.txt	--	-

ソースコード

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 shuffle
    p_1_pfs = prime_factors(p-1, return_e=False)
    candidates_for_r = list(range(2, p))
    shuffle(candidates_for_r)  # 乱択で原始根を探す
    for cfr in 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()