結果

問題 No.981 一般冪乗根
ユーザー hari64hari64
提出日時 2021-07-22 15:07:35
言語 Python3
(3.12.2 + numpy 1.26.4 + scipy 1.12.0)
結果
RE  
実行時間 -
コード長 8,272 bytes
コンパイル時間 82 ms
コンパイル使用メモリ 13,184 KB
実行使用メモリ 65,876 KB
最終ジャッジ日時 2024-07-17 14:47:26
合計ジャッジ時間 12,980 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge2
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
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 -- -
権限があれば一括ダウンロードができます

ソースコード

diff # 

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  # 原始根でない限り、必ず存在するとは言えない


"""
関連する話題 これはメモ20210505

フェルマーの小定理
pが素数 aがpの倍数でない正の整数の時 a^(p-1)≡1 (mod p)
証明 a^p=aを示す 帰納法にすると(1+a)^p≡1+a^p≡1+a (∵二項係数) よって示された
別解 ia(i={1,2,…,p-1})にはmod pで1,2,…,p-1が出る(∵背理法)
     よってa^(p-1)*(p-1)!≡(p-1)!⇔a^(p-1)≡1(∵(p-1)!とpは互いに素)
https://manabitimes.jp/math/680
https://qiita.com/drken/items/6b4031ccbb2cab7436f3

レプユニット数(11…11) これは完全に数学 https://manabitimes.jp/math/1374

NTT(整数環FFT)(Number Theoretic Transform)
FFTで出てくる丸め誤差の問題を解決する畳み込みの手法
https://math314.hateblo.jp/entry/2015/05/07/014908

garnerのアルゴリズム 本質的にはあまり変わらない
https://www.csee.umbc.edu/~lomonaco/s08/441/handouts/Garner-Alg-Example.pdf
"""


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 crt(rs: list, mods: list) -> tuple:
    """
    同じ長さのリストr,mを引数に取り、このリストの長さをnとした時
    x≡rs[i] (mod ms[i]) ∀i∊{0,1,…,n-1}を解く
    答えが存在するならばx≡y(mod z) (0<=y<z=lcm(ms))なる(y,z)を返す
    制約 len(rs)=len(ms) 1<=ms[i] 法が互いに素である必要はない
    計算量 O(nlog(lcm(ms)))
    簡単に言えば中国剰余定理 (Chinese Remainder Theorem)のことである
    https://qiita.com/drken/items/ae02240cd1f8edfc86fd
    https://github.com/atcoder/ac-library/blob/master/atcoder/math.hpp
    https://manabitimes.jp/math/837 https://manabitimes.jp/math/838
    中国剰余定理は、一次不定方程式(ベズー等式)を繰り返し解くことにほぼ同義
    ax+by=cが整数解を持つ必要十分条件がcがgcd(a,b)の倍数であることなどを思い出しておくと分かりやすい
    https://manabitimes.jp/math/674
    """
    assert len(rs) == len(mods)
    n = len(rs)
    r0 = 0  # 0<=r0<m0という条件の下で考えていく
    m0 = 1  # x≡0 (mod 1)は恒等式 これにより最初のステップを例外扱いしないで済む
    for i in range(n):  # x≡r0(mod m0)≡r1(mod m1)という連立式を繰り返し解く
        assert 1 <= mods[i]
        r1 = rs[i] % mods[i]
        m1 = mods[i]
        if m0 < m1:  # m1<=m0にする
            r0, r1 = r1, r0
            m0, m1 = m1, m0
        if m0 % m1 == 0:
            if r0 % m1 != r1:  # 大小関係より矛盾
                return (0, 0)
            else:  # 新しい式の主張は既に満たされている
                continue
        # 以上の処理でm1<m0,2*max(m0,m1)<=lcm(m0,m1)が保証される
        # 解をx≡r2 (mod lcm(m0,m1))と置くと
        # r2≡r0 (mod m0) r2≡r1 (mod m1) より (r0+x*m0)≡r1 (mod m1)
        # つまり x*m0≡r1-r0 (mod m1) これを法も含めてg(=gcd(m0,m1))で割ると
        # x≡((r1-r0)*(u0^-1))/g (mod u1)  (m0/g=u0 m1/g=u1)
        # この辺りは少しACLと説明を変えている そちらもそちらで分かりやすい
        g, im = _inv_gcd(m0, m1)  # im=(u0)^-1 (mod u1) (0<=im<u1)
        u1 = m1 // g
        if (r1 - r0) % g:  # 互いにその場合に帰着できないならそれは解けない
            return (0, 0)
        x = (r1 - r0) // g * im % u1  # 先述の式
        r0 += x * m0  # ≡r2 (mod m2)
        m0 *= u1  # =m2(=lcm(m0,m1))
        if r0 < 0:
            r0 += m0
    return (r0, m0)


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