half-GCD
code:halfgcd.sage.py
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
"""
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)
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)
r = b0 % b1
if r == 0:
return b1
return pgcd(b1, r)
Pari/GP にはhalf gcdが実装されているのでそれを使うこともできる code:sage
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')))
参考文献
を読もうー完ー