がたがた言わずに与えられたものを使え。これは 一変数多項式のgcd の高速なやつです
def pdivmod(u, v): """ polynomial version of divmod """ q = u // v r = u - q*v return (q, r) def hgcd(u, v, min_degree=10): """ Calculate Half-GCD of (u, v) f and g are univariate polynomial http://web.cs.iastate.edu/~cs577/handouts/polydivide.pdf """ x = u.parent().gen() if u.degree() < v.degree(): u, v = v, u if 2*v.degree() < u.degree() or u.degree() < min_degree: q = u // v return matrix([[1, -q], [0, 1]]) m = u.degree() // 2 b0, c0 = pdivmod(u, x^m) b1, c1 = pdivmod(v, x^m) R = hgcd(b0, b1) DE = R * matrix([[u], [v]]) d, e = DE[0,0], DE[1,0] q, f = pdivmod(d, e) g0 = e // x^(m//2) g1 = f // x^(m//2) S = hgcd(g0, g1) return S * matrix([[0, 1], [1, -q]]) * R def pgcd(u, v): """ fast implementation of polynomial GCD using hgcd """ if u.degree() < v.degree(): u, v = v, u if v == 0: return u if u % v == 0: return v if u.degree() < 10: while v != 0: u, v = v, u % v return u R = hgcd(u, v) B = R * matrix([[u], [v]]) b0, b1 = B[0,0], B[1,0] r = b0 % b1 if r == 0: return b1 return pgcd(b1, r)
Pari/GP にはhalf gcdが実装されているのでそれを使うこともできる
PR.<x> = PolynomialRing(Zmod(N)) f1 = x**e - c1 f2 = (x + a)**e - c2 g = PR(f1._pari_with_name('x').gcd(f2._pari_with_name('x')))
よすぽジャッジのinv of polynomials ( https://judge.yosupo.jp/problem/inv_of_polynomials ) の解説。
参考文献
(1) https://dl.acm.org/doi/10.1145/800125.804045
(2) https://pdfs.semanticscholar.org/a7e7/b01a3dd6ac0ec160b35e513c5efa38c2369e.pdf
を読もうー完ー