half-GCD

がたがた言わずに与えられたものを使え。これは 一変数多項式の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

を読もうー完ー

https://scrapbox.io/37zigen-43465887/half-GCD