結果
| 問題 |
No.981 一般冪乗根
|
| ユーザー |
 hari64
|
| 提出日時 | 2021-07-22 15:02:58 |
| 言語 | PyPy3 (7.3.15) |
| 結果 |
RE
|
| 実行時間 | - |
| コード長 | 8,272 bytes |
| コンパイル時間 | 150 ms |
| コンパイル使用メモリ | 82,032 KB |
| 実行使用メモリ | 271,736 KB |
| 最終ジャッジ日時 | 2024-07-17 14:46:41 |
| 合計ジャッジ時間 | 13,915 ms |
|
ジャッジサーバーID （参考情報） |
judge5 / judge2 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | RE * 5 TLE * 1 -- * 38 |
ソースコード
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()