結果

問題 No.981 一般冪乗根
ユーザー gew1fw
提出日時 2025-06-12 13:28:50
言語 PyPy3
(7.3.15)
結果
WA  
実行時間 -
コード長 4,204 bytes
コンパイル時間 414 ms
コンパイル使用メモリ 82,216 KB
実行使用メモリ 77,232 KB
最終ジャッジ日時 2025-06-12 13:34:35
合計ジャッジ時間 41,918 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge4
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 6 WA * 38
権限があれば一括ダウンロードができます

ソースコード

diff # 

import sys
import random
from math import gcd

def input():
    return sys.stdin.read()

def extended_gcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = extended_gcd(b % a, a)
        return (g, x - (b // a) * y, y)

def inverse(a, m):
    g, x, y = extended_gcd(a, m)
    if g != 1:
        return None
    else:
        return x % m

def tonelli_shanks(n, p):
    if n == 0:
        return 0
    ls = pow(n, (p-1)//2, p)
    if ls != 1:
        return -1
    if p % 4 == 3:
        x = pow(n, (p+1)//4, p)
        return x
    Q = p-1
    S = 0
    while Q % 2 == 0:
        Q //= 2
        S += 1
    z = 2
    while pow(z, (p-1)//2, p) != p-1:
        z += 1
    c = pow(z, Q, p)
    x = pow(n, (Q+1)//2, p)
    t = pow(n, Q, p)
    m = S
    while t != 1:
        i, temp = 0, t
        while temp != 1 and i < m:
            temp = pow(temp, 2, p)
            i += 1
        if i == m:
            return -1
        b = pow(c, 2**(m - i -1), p)
        x = (x * b) % p
        t = (t * b * b) % p
        c = (b * b) % p
        m = i
    return x

def discrete_log_pollard(h, delta, p, ord):
    def new_xab(x, a, b):
        subset = x % 3
        if subset == 0:
            x = (x * x) % p
            a = (a * 2) % ord
            b = (b * 2) % ord
        elif subset == 1:
            x = (x * h) % p
            a = (a + 1) % ord
        else:
            x = (x * delta) % p
            b = (b + 1) % ord
        return x, a, b

    for _ in range(100):
        x, a, b = 1, 0, 0
        X, A, B = x, a, b
        for _ in range(1, ord):
            x, a, b = new_xab(x, a, b)
            X, A, B = new_xab(X, A, B)
            X, A, B = new_xab(X, A, B)
            if x == X:
                break
        else:
            continue
        denominator = (a - A) % ord
        numerator = (B - b) % ord
        if denominator == 0:
            if numerator == 0:
                return 0
            else:
                return -1
        g, x_gcd, y_gcd = extended_gcd(denominator, ord)
        if numerator % g != 0:
            return -1
        x0 = ((numerator // g) * x_gcd) % (ord // g)
        for i in range(g):
            m_candidate = (x0 + i * (ord // g)) % ord
            if pow(h, m_candidate, p) == delta:
                return m_candidate
        return -1
    return -1

def amm_root(a, e, p):
    if a == 0:
        return 0
    t = 0
    s = p-1
    while s % e == 0:
        s //= e
        t += 1
    if t == 0:
        inv_e = inverse(e, s)
        return pow(a, inv_e, p) if inv_e else -1
    g = 2
    while pow(g, s, p) == 1:
        g += 1
    h = pow(g, s, p)
    inv_e_t = inverse(e**t, s)
    if not inv_e_t:
        return -1
    a0 = pow(a, inv_e_t, p)
    gamma = e ** (t-1)
    inv_e = inverse(e, gamma)
    if not inv_e:
        return -1
    x0 = pow(a0, inv_e, p)
    c = pow(g, s * gamma, p)
    h_power = pow(h, gamma, p)
    current_e = e
    for i in range(1, t):
        delta = (a0 * inverse(pow(x0, e, p), p)) % p
        delta = pow(delta, current_e // e, p)
        current_e = current_e // e
        m = discrete_log_pollard(h_power, delta, p, e)
        if m == -1:
            return -1
        x0 = (x0 * pow(c, m, p)) % p
        c = pow(c, e, p)
        h_power = pow(h_power, e, p)
    return x0

def solve():
    data = input().split()
    T = int(data[0])
    idx = 1
    for _ in range(T):
        p = int(data[idx])
        k = int(data[idx+1])
        a = int(data[idx+2])
        idx +=3
        if a == 0:
            print(0)
            continue
        d = gcd(k, p-1)
        exponent = (p-1) // d
        if pow(a, exponent, p) != 1:
            print(-1)
            continue
        if d == 1:
            y = a
        elif d == 2:
            y = tonelli_shanks(a, p)
            if y == -1:
                print(-1)
                continue
        else:
            y = amm_root(a, d, p)
            if y == -1:
                print(-1)
                continue
        k_prime = k // d
        m = (p-1) // d
        inv_k_prime = inverse(k_prime, m)
        if not inv_k_prime:
            print(-1)
            continue
        x = pow(y, inv_k_prime, p)
        print(x)

solve()
0