# Source code for diofant.polys.euclidtools

"""Euclidean algorithms, GCDs, LCMs and polynomial remainder sequences. """

from ..core import cacheit
from ..ntheory import nextprime
from ..ntheory.modular import crt, symmetric_residue
from .densearith import (dmp_add, dmp_div, dmp_max_norm, dmp_mul,
dmp_mul_ground, dmp_mul_term, dmp_neg, dmp_pow,
dmp_quo, dmp_quo_ground, dmp_rem, dmp_sub,
dmp_sub_mul, dup_mul)
from .densebasic import (dmp_apply_pairs, dmp_convert, dmp_deflate,
dmp_degree_in, dmp_ground, dmp_ground_LC, dmp_inflate,
dmp_LC, dmp_one, dmp_one_p, dmp_raise, dmp_strip,
dmp_zero, dmp_zero_p, dmp_zeros)
from .densetools import (dmp_clear_denoms, dmp_eval_in, dmp_ground_monic,
dmp_ground_primitive, dmp_ground_trunc)
from .heuristicgcd import heugcd
from .polyconfig import query
from .polyerrors import (DomainError, HeuristicGCDFailed, HomomorphismFailed,
NotInvertible)

[docs]def dup_half_gcdex(f, g, K):
"""
Half extended Euclidean algorithm in F[x].

Returns (s, h) such that h = gcd(f, g) and s*f = h (mod g).

Examples
========

>>> R, x = ring("x", QQ)

>>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15
>>> g = x**3 + x**2 - 4*x - 4

>>> f.half_gcdex(g)
(-1/5*x + 3/5, x + 1)

"""
if not K.is_Field:
raise DomainError("can't compute half extended GCD over %s" % K)

a, b = [K.one], []

while g:
q, r = dmp_div(f, g, 0, K)
f, g = g, r
a, b = b, dmp_sub_mul(a, q, b, 0, K)

a = dmp_quo_ground(a, dmp_LC(f, K), 0, K)
f = dmp_ground_monic(f, 0, K)

return a, f

[docs]def dup_gcdex(f, g, K):
"""
Extended Euclidean algorithm in F[x].

Returns (s, t, h) such that h = gcd(f, g) and s*f + t*g = h.

Examples
========

>>> R, x = ring("x", QQ)

>>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15
>>> g = x**3 + x**2 - 4*x - 4

>>> f.gcdex(g)
(-1/5*x + 3/5, 1/5*x**2 - 6/5*x + 2, x + 1)

"""
s, h = dup_half_gcdex(f, g, K)

F = dmp_sub_mul(h, s, f, 0, K)
t = dmp_quo(F, g, 0, K)

return s, t, h

def dup_invert(f, g, K):
"""
Compute multiplicative inverse of f modulo g in F[x].

Examples
========

>>> R, x = ring("x", QQ)

>>> f = x**2 - 1
>>> g = 2*x - 1
>>> h = x - 1

>>> R.dup_invert(f, g)
-4/3

>>> R.dup_invert(f, h)
Traceback (most recent call last):
...
NotInvertible: zero divisor

"""
s, h = dup_half_gcdex(f, g, K)

if h == [K.one]:
return dmp_rem(s, g, 0, K)
else:
raise NotInvertible("zero divisor")

def dmp_prem(f, g, u, K):
"""
Polynomial pseudo-remainder in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> (x**2 + x*y).prem(2*x + 2)
-4*y + 4

References
==========

* :cite:Knuth1985seminumerical, p. 407.

"""
df = dmp_degree_in(f, 0, u)
dg = dmp_degree_in(g, 0, u)

if dg < 0:
raise ZeroDivisionError("polynomial division")

r, dr = f, df

if df < dg:
return r

N = df - dg + 1
lc_g = dmp_LC(g, K)

while True:
lc_r = dmp_LC(r, K)
j, N = dr - dg, N - 1

R = dmp_mul_term(r, lc_g, 0, u, K)
G = dmp_mul_term(g, lc_r, j, u, K)
r = dmp_sub(R, G, u, K)

_dr, dr = dr, dmp_degree_in(r, 0, u)

if dr < dg:
break
assert dr < _dr

if u:
c = dmp_pow(lc_g, N, u - 1, K)
else:
c = lc_g**N

return dmp_mul_term(r, c, 0, u, K)

def dup_inner_subresultants(f, g, K):
"""
Subresultant PRS algorithm in K[x].

Computes the subresultant polynomial remainder sequence (PRS)
and the non-zero scalar subresultants of f and g.
By [1] Thm. 3, these are the constants '-c' (- to optimize
computation of sign).
The first subdeterminant is set to 1 by convention to match
the polynomial and the scalar subdeterminants.
If 'deg(f) < deg(g)', the subresultants of '(g,f)' are computed.

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dmp_inner_subresultants(x**2 + 1, x**2 - 1)
([x**2 + 1, x**2 - 1, -2], [1, 1, 4])

References
==========

* :cite:Brown1978prs

"""
n = dmp_degree_in(f, 0, 0)
m = dmp_degree_in(g, 0, 0)

if n < m:
f, g = g, f
n, m = m, n

if not f:
return [], []

if not g:
return [f], [K.one]

R = [f, g]
d = n - m

b = (-K.one)**(d + 1)

h = dmp_prem(f, g, 0, K)
h = dmp_mul_ground(h, b, 0, K)

lc = dmp_LC(g, K)
c = lc**d

# Conventional first scalar subdeterminant is 1
S = [K.one, c]
c = -c

while h:
k = dmp_degree_in(h, 0, 0)
R.append(h)

f, g, m, d = g, h, k, m - k

b = -lc * c**d

h = dmp_prem(f, g, 0, K)
h = dmp_quo_ground(h, b, 0, K)

lc = dmp_LC(g, K)

if d > 1:        # abnormal case
q = c**(d - 1)
c = K.quo((-lc)**d, q)
else:
c = -lc

S.append(-c)

return R, S

def dup_prs_resultant(f, g, K):
"""
Resultant algorithm in K[x] using subresultant PRS.

Examples
========

>>> R, x = ring("x", ZZ)

>>> (x**2 + 1).resultant(x**2 - 1, includePRS=True)
(4, [x**2 + 1, x**2 - 1, -2])

"""
if not f or not g:
return K.zero, []

R, S = dup_inner_subresultants(f, g, K)

if dmp_degree_in(R[-1], 0, 0) > 0:
return K.zero, R

return S[-1], R

[docs]def dmp_inner_subresultants(f, g, u, K):
"""
Subresultant PRS algorithm in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = 3*x**2*y - y**3 - 4
>>> g = x**2 + x*y**3 - 9

>>> a = 3*x*y**4 + y**3 - 27*y + 4
>>> b = -3*y**10 - 12*y**7 + y**6 - 54*y**4 + 8*y**3 + 729*y**2 - 216*y + 16

>>> prs = [f, g, a, b]
>>> sres = [[1], [1], [3, 0, 0, 0, 0], [-3, 0, 0, -12, 1, 0, -54, 8, 729, -216, 16]]

>>> R.dmp_inner_subresultants(f, g) == (prs, sres)
True

"""
if not u:
return dup_inner_subresultants(f, g, K)

n = dmp_degree_in(f, 0, u)
m = dmp_degree_in(g, 0, u)

if n < m:
f, g = g, f
n, m = m, n

if dmp_zero_p(f, u):
return [], []

v = u - 1
if dmp_zero_p(g, u):
return [f], [dmp_ground(K.one, v)]

R = [f, g]
d = n - m

b = dmp_pow(dmp_ground(-K.one, v), d + 1, v, K)

h = dmp_prem(f, g, u, K)
h = dmp_mul_term(h, b, 0, u, K)

lc = dmp_LC(g, K)
c = dmp_pow(lc, d, v, K)

S = [dmp_ground(K.one, v), c]
c = dmp_neg(c, v, K)

while not dmp_zero_p(h, u):
k = dmp_degree_in(h, 0, u)
R.append(h)

f, g, m, d = g, h, k, m - k

b = dmp_mul(dmp_neg(lc, v, K),
dmp_pow(c, d, v, K), v, K)

h = dmp_prem(f, g, u, K)
h = [dmp_quo(ch, b, v, K) for ch in h]

lc = dmp_LC(g, K)

if d > 1:
p = dmp_pow(dmp_neg(lc, v, K), d, v, K)
q = dmp_pow(c, d - 1, v, K)
c = dmp_quo(p, q, v, K)
else:
c = dmp_neg(lc, v, K)

S.append(dmp_neg(c, v, K))

return R, S

[docs]def dmp_subresultants(f, g, u, K):
"""
Computes subresultant PRS of two polynomials in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = 3*x**2*y - y**3 - 4
>>> g = x**2 + x*y**3 - 9

>>> a = 3*x*y**4 + y**3 - 27*y + 4
>>> b = -3*y**10 - 12*y**7 + y**6 - 54*y**4 + 8*y**3 + 729*y**2 - 216*y + 16

>>> f.subresultants(g) == [f, g, a, b]
True

"""
return dmp_inner_subresultants(f, g, u, K)[0]

def dmp_prs_resultant(f, g, u, K):
"""
Resultant algorithm in K[X] using subresultant PRS.

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = 3*x**2*y - y**3 - 4
>>> g = x**2 + x*y**3 - 9

>>> a = 3*x*y**4 + y**3 - 27*y + 4
>>> b = -3*y**10 - 12*y**7 + y**6 - 54*y**4 + 8*y**3 + 729*y**2 - 216*y + 16

>>> res, prs = R.dmp_resultant(f, g, includePRS=True)

>>> res == b             # resultant has n-1 variables
False
>>> res == b.drop(x)
True
>>> prs == [f, g, a, b]
True

"""
if not u:
return dup_prs_resultant(f, g, K)

if dmp_zero_p(f, u) or dmp_zero_p(g, u):
return dmp_zero(u - 1), []

R, S = dmp_inner_subresultants(f, g, u, K)

if dmp_degree_in(R[-1], 0, u) > 0:
return dmp_zero(u - 1), R

return S[-1], R

[docs]def dmp_zz_modular_resultant(f, g, p, u, K):
"""
Compute resultant of f and g modulo a prime p.

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x + y + 2
>>> g = 2*x*y + x + 3

>>> R.dmp_zz_modular_resultant(f, g, 5)
-2*y**2 + 1

"""
if not u:
return symmetric_residue(dup_prs_resultant(f, g, K)[0] % p, p)

v = u - 1

n = dmp_degree_in(f, 0, u)
m = dmp_degree_in(g, 0, u)

N = dmp_degree_in(f, 1, u)
M = dmp_degree_in(g, 1, u)

B = n*M + m*N

D, a = [K.one], -K.one
r = dmp_zero(v)

while dmp_degree_in(D, 0, 0) <= B:
while True:
a += K.one

if a == p:
raise HomomorphismFailed('no luck')

F = dmp_eval_in(f, symmetric_residue(a, p), 1, u, K)

if dmp_degree_in(F, 0, v) == n:
G = dmp_eval_in(g, symmetric_residue(a, p), 1, u, K)

if dmp_degree_in(G, 0, v) == m:
break

R = dmp_zz_modular_resultant(F, G, p, v, K)
e = dmp_eval_in(r, a, 0, v, K)

if not v:
R = dmp_strip([R], 0)
e = dmp_strip([e], 0)
else:
R = [R]
e = [e]

d = K.invert(dmp_eval_in(D, a, 0, 0, K), p)
d = dmp_mul_ground(D, d, 0, K)
d = dmp_raise(d, v, 0, K)

c = dmp_mul(d, dmp_sub(R, e, v, K), v, K)
r = dmp_add(r, c, v, K)

r = dmp_ground_trunc(r, p, v, K)

D = dup_mul(D, [K.one, -a], K)
D = dmp_ground_trunc(D, p, 0, K)

return r

def _collins_crt(r, R, P, p, K):
"""Wrapper of CRT for Collins's resultant algorithm."""
return K(crt([P, p], [r, R], check=False, symmetric=True)[0])

[docs]def dmp_zz_collins_resultant(f, g, u, K):
"""
Collins's modular resultant algorithm in Z[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x + y + 2
>>> g = 2*x*y + x + 3

>>> R.dmp_zz_collins_resultant(f, g)
-2*y**2 - 5*y + 1

"""

n = dmp_degree_in(f, 0, u)
m = dmp_degree_in(g, 0, u)

if n < 0 or m < 0:
return dmp_zero(u - 1)

A = dmp_max_norm(f, u, K)
B = dmp_max_norm(g, u, K)

a = dmp_ground_LC(f, u, K)
b = dmp_ground_LC(g, u, K)

v = u - 1

B = K(2)*K.factorial(K(n + m))*A**m*B**n
r, p, P = dmp_zero(v), K.one, K.one

while P <= B:
p = K(nextprime(p))

while not (a % p) or not (b % p):
p = K(nextprime(p))

F = dmp_ground_trunc(f, p, u, K)
G = dmp_ground_trunc(g, p, u, K)

try:
R = dmp_zz_modular_resultant(F, G, p, u, K)
except HomomorphismFailed:
continue

if P == K.one:
r = R
else:
r = dmp_apply_pairs(r, R, _collins_crt, (P, p, K), v, K)

P *= p

return r

[docs]def dmp_qq_collins_resultant(f, g, u, K0):
"""
Collins's modular resultant algorithm in Q[X].

Examples
========

>>> R, x, y = ring("x y", QQ)

>>> f = x/2 + y + QQ(2, 3)
>>> g = 2*x*y + x + 3

>>> R.dmp_qq_collins_resultant(f, g)
-2*y**2 - 7/3*y + 5/6

"""
n = dmp_degree_in(f, 0, u)
m = dmp_degree_in(g, 0, u)

if n < 0 or m < 0:
return dmp_zero(u - 1)

K1 = K0.ring

cf, f = dmp_clear_denoms(f, u, K0, K1)
cg, g = dmp_clear_denoms(g, u, K0, K1)

f = dmp_convert(f, u, K0, K1)
g = dmp_convert(g, u, K0, K1)

r = dmp_zz_collins_resultant(f, g, u, K1)
r = dmp_convert(r, u - 1, K1, K0)

c = K0.convert(cf**m * cg**n, K1)

return dmp_quo_ground(r, c, u - 1, K0)

[docs]@cacheit
def dmp_resultant(f, g, u, K, includePRS=False):
"""
Computes resultant of two polynomials in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = 3*x**2*y - y**3 - 4
>>> g = x**2 + x*y**3 - 9

>>> f.resultant(g)
-3*y**10 - 12*y**7 + y**6 - 54*y**4 + 8*y**3 + 729*y**2 - 216*y + 16

"""
if includePRS:
return dmp_prs_resultant(f, g, u, K)

if K.is_Field:
if K.is_RationalField and query('USE_COLLINS_RESULTANT'):
return dmp_qq_collins_resultant(f, g, u, K)
else:
if K.is_IntegerRing and query('USE_COLLINS_RESULTANT'):
return dmp_zz_collins_resultant(f, g, u, K)

return dmp_prs_resultant(f, g, u, K)[0]

def _dmp_rr_trivial_gcd(f, g, u, K):
"""Handle trivial cases in GCD algorithm over a ring."""
zero_f = dmp_zero_p(f, u)
zero_g = dmp_zero_p(g, u)

if zero_f and zero_g:
return tuple(dmp_zeros(3, u, K))
elif zero_f:
if K.is_negative(dmp_ground_LC(g, u, K)):
return dmp_neg(g, u, K), dmp_zero(u), dmp_ground(-K.one, u)
else:
return g, dmp_zero(u), dmp_one(u, K)
elif zero_g:
if K.is_negative(dmp_ground_LC(f, u, K)):
return dmp_neg(f, u, K), dmp_ground(-K.one, u), dmp_zero(u)
else:
return f, dmp_one(u, K), dmp_zero(u)
elif dmp_one_p(f, u, K) or dmp_one_p(g, u, K):
return dmp_one(u, K), f, g
elif u and query('USE_SIMPLIFY_GCD'):
return _dmp_simplify_gcd(f, g, u, K)

def _dmp_ff_trivial_gcd(f, g, u, K):
"""Handle trivial cases in GCD algorithm over a field."""
zero_f = dmp_zero_p(f, u)
zero_g = dmp_zero_p(g, u)

if zero_f and zero_g:
return tuple(dmp_zeros(3, u, K))
elif zero_f:
return (dmp_ground_monic(g, u, K),
dmp_zero(u),
dmp_ground(dmp_ground_LC(g, u, K), u))
elif zero_g:
return (dmp_ground_monic(f, u, K),
dmp_ground(dmp_ground_LC(f, u, K), u),
dmp_zero(u))
elif u and query('USE_SIMPLIFY_GCD'):
return _dmp_simplify_gcd(f, g, u, K)

def _dmp_simplify_gcd(f, g, u, K):
"""Try to eliminate x_0 from GCD computation in K[X]."""
df = dmp_degree_in(f, 0, u)
dg = dmp_degree_in(g, 0, u)

if df > 0 and dg > 0:
return

if not (df or dg):
F = dmp_LC(f, K)
G = dmp_LC(g, K)
else:
if not df:
F = dmp_LC(f, K)
G = dmp_content(g, u, K)
else:
F = dmp_content(f, u, K)
G = dmp_LC(g, K)

v = u - 1
h = dmp_gcd(F, G, v, K)

cff = [dmp_quo(cf, h, v, K) for cf in f]
cfg = [dmp_quo(cg, h, v, K) for cg in g]

return [h], cff, cfg

def dup_rr_prs_gcd(f, g, K):
"""
Computes polynomial GCD using subresultants over a ring.

Returns (h, cff, cfg) such that a = gcd(f, g), cff = quo(f, h),
and cfg = quo(g, h).

Examples
========

>>> R, x = ring("x", ZZ)

>>> R.dup_rr_prs_gcd(x**2 - 1, x**2 - 3*x + 2)
(x - 1, x + 1, x - 2)

"""
result = _dmp_rr_trivial_gcd(f, g, 0, K)

if result is not None:
return result

fc, F = dmp_ground_primitive(f, 0, K)
gc, G = dmp_ground_primitive(g, 0, K)

c = K.gcd(fc, gc)

h = dmp_subresultants(F, G, 0, K)[-1]
_, h = dmp_ground_primitive(h, 0, K)

h = dmp_mul_ground(h, c, 0, K)

cff = dmp_quo(f, h, 0, K)
cfg = dmp_quo(g, h, 0, K)

return h, cff, cfg

def dup_ff_prs_gcd(f, g, K):
"""
Computes polynomial GCD using subresultants over a field.

Returns (h, cff, cfg) such that a = gcd(f, g), cff = quo(f, h),
and cfg = quo(g, h).

Examples
========

>>> R, x = ring("x", QQ)

>>> R.dup_ff_prs_gcd(x**2 - 1, x**2 - 3*x + 2)
(x - 1, x + 1, x - 2)

"""
result = _dmp_ff_trivial_gcd(f, g, 0, K)

if result is not None:
return result

h = dmp_subresultants(f, g, 0, K)[-1]
h = dmp_ground_monic(h, 0, K)

cff = dmp_quo(f, h, 0, K)
cfg = dmp_quo(g, h, 0, K)

return h, cff, cfg

[docs]def dmp_rr_prs_gcd(f, g, u, K):
"""
Computes polynomial GCD using subresultants over a ring.

Returns (h, cff, cfg) such that a = gcd(f, g), cff = quo(f, h),
and cfg = quo(g, h).

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x**2 + 2*x*y + y**2
>>> g = x**2 + x*y

>>> R.dmp_rr_prs_gcd(f, g)
(x + y, x + y, x)

"""
if not u:
return dup_rr_prs_gcd(f, g, K)

result = _dmp_rr_trivial_gcd(f, g, u, K)

if result is not None:
return result

fc, F = dmp_primitive(f, u, K)
gc, G = dmp_primitive(g, u, K)

h = dmp_subresultants(F, G, u, K)[-1]
c, _, _ = dmp_rr_prs_gcd(fc, gc, u - 1, K)

if K.is_negative(dmp_ground_LC(h, u, K)):
h = dmp_neg(h, u, K)

_, h = dmp_primitive(h, u, K)
h = dmp_mul_term(h, c, 0, u, K)

cff = dmp_quo(f, h, u, K)
cfg = dmp_quo(g, h, u, K)

return h, cff, cfg

[docs]def dmp_ff_prs_gcd(f, g, u, K):
"""
Computes polynomial GCD using subresultants over a field.

Returns (h, cff, cfg) such that a = gcd(f, g), cff = quo(f, h),
and cfg = quo(g, h).

Examples
========

>>> R, x, y = ring("x y", QQ)

>>> f = x**2/2 + x*y + y**2/2
>>> g = x**2 + x*y

>>> R.dmp_ff_prs_gcd(f, g)
(x + y, 1/2*x + 1/2*y, x)

"""
if not u:
return dup_ff_prs_gcd(f, g, K)

result = _dmp_ff_trivial_gcd(f, g, u, K)

if result is not None:
return result

fc, F = dmp_primitive(f, u, K)
gc, G = dmp_primitive(g, u, K)

h = dmp_subresultants(F, G, u, K)[-1]
c, _, _ = dmp_ff_prs_gcd(fc, gc, u - 1, K)

_, h = dmp_primitive(h, u, K)
h = dmp_mul_term(h, c, 0, u, K)
h = dmp_ground_monic(h, u, K)

cff = dmp_quo(f, h, u, K)
cfg = dmp_quo(g, h, u, K)

return h, cff, cfg

[docs]def dmp_zz_heu_gcd(f, g, u, K):
"""
Heuristic polynomial GCD in Z[X].

Returns (h, cff, cfg) such that a = gcd(f, g),
cff = quo(f, h), and cfg = quo(g, h).

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x**2 + 2*x*y + y**2
>>> g = x**2 + x*y

>>> R.dmp_zz_heu_gcd(f, g)
(x + y, x + y, x)

See Also
========

diofant.polys.heuristicgcd.heugcd

References
==========

* :cite:Liao1995heuristic

"""
result = _dmp_rr_trivial_gcd(f, g, u, K)

if result is not None:
return result

ring = K.poly_ring(*["_%d" % i for i in range(u + 1)])
f, g = map(ring.from_dense, (f, g))
return tuple(map(ring.to_dense, heugcd(f, g)))

[docs]def dmp_qq_heu_gcd(f, g, u, K0):
"""
Heuristic polynomial GCD in Q[X].

Returns (h, cff, cfg) such that a = gcd(f, g),
cff = quo(f, h), and cfg = quo(g, h).

Examples
========

>>> R, x, y = ring("x y", QQ)

>>> f = x**2/4 + x*y + y**2
>>> g = x**2/2 + x*y

>>> R.dmp_qq_heu_gcd(f, g)
(x + 2*y, 1/4*x + 1/2*y, 1/2*x)

"""
result = _dmp_ff_trivial_gcd(f, g, u, K0)

if result is not None:
return result

K1 = K0.ring

cf, f = dmp_clear_denoms(f, u, K0, K1)
cg, g = dmp_clear_denoms(g, u, K0, K1)

f = dmp_convert(f, u, K0, K1)
g = dmp_convert(g, u, K0, K1)

h, cff, cfg = dmp_zz_heu_gcd(f, g, u, K1)

h = dmp_convert(h, u, K1, K0)

c = dmp_ground_LC(h, u, K0)
h = dmp_ground_monic(h, u, K0)

cff = dmp_convert(cff, u, K1, K0)
cfg = dmp_convert(cfg, u, K1, K0)

cff = dmp_mul_ground(cff, K0.quo(c, cf), u, K0)
cfg = dmp_mul_ground(cfg, K0.quo(c, cg), u, K0)

return h, cff, cfg

def _dmp_zz_modgcd(f, g, u, K):
from .modulargcd import modgcd
ring = K.poly_ring(*["_%d" % i for i in range(u + 1)])
f, g = map(ring.from_dense, (f, g))
return tuple(map(ring.to_dense, modgcd(f, g)))

_gcd_zz_methods = {'modgcd': _dmp_zz_modgcd,
'prs': dmp_rr_prs_gcd}

def _dmp_aa_modgcd(f, g, u, K):
from .modulargcd import func_field_modgcd
ring = K.poly_ring(*["_%d" % i for i in range(u + 1)])
f, g = map(ring.from_dense, (f, g))
return tuple(map(ring.to_dense, func_field_modgcd(f, g)))

_gcd_aa_methods = {'modgcd': _dmp_aa_modgcd,
'prs': dmp_ff_prs_gcd}

def _dmp_inner_gcd(f, g, u, K):
"""Helper function for dmp_inner_gcd()."""
if not K.is_Exact:
try:
exact = K.get_exact()
except DomainError:
return dmp_one(u, K), f, g

f = dmp_convert(f, u, K, exact)
g = dmp_convert(g, u, K, exact)

h, cff, cfg = _dmp_inner_gcd(f, g, u, exact)

h = dmp_convert(h, u, exact, K)
cff = dmp_convert(cff, u, exact, K)
cfg = dmp_convert(cfg, u, exact, K)

return h, cff, cfg
elif K.is_Field:
if K.is_RationalField:
if query('USE_HEU_GCD'):
try:
return dmp_qq_heu_gcd(f, g, u, K)
except HeuristicGCDFailed:  # pragma: no cover
pass
elif K.is_AlgebraicField:
method = _gcd_aa_methods[query('GCD_AA_METHOD')]
return method(f, g, u, K)

return dmp_ff_prs_gcd(f, g, u, K)
else:
if K.is_IntegerRing:
if query('USE_HEU_GCD'):
try:
return dmp_zz_heu_gcd(f, g, u, K)
except HeuristicGCDFailed:  # pragma: no cover
pass

method = _gcd_zz_methods[query('FALLBACK_GCD_ZZ_METHOD')]
return method(f, g, u, K)

return dmp_rr_prs_gcd(f, g, u, K)

[docs]def dmp_inner_gcd(f, g, u, K):
"""
Computes polynomial GCD and cofactors of f and g in K[X].

Returns (h, cff, cfg) such that h = gcd(f, g),
cff = quo(f, h), and cfg = quo(g, h).

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x**2 + 2*x*y + y**2
>>> g = x**2 + x*y

>>> R.dmp_inner_gcd(f, g)
(x + y, x + y, x)

"""
J, (f, g) = dmp_deflate((f, g), u, K)
h, cff, cfg = _dmp_inner_gcd(f, g, u, K)

return (dmp_inflate(h, J, u, K),
dmp_inflate(cff, J, u, K),
dmp_inflate(cfg, J, u, K))

[docs]def dmp_gcd(f, g, u, K):
"""
Computes polynomial GCD of f and g in K[X].

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> f = x**2 + 2*x*y + y**2
>>> g = x**2 + x*y

>>> R.dmp_gcd(f, g)
x + y

"""
return dmp_inner_gcd(f, g, u, K)[0]

[docs]def dmp_content(f, u, K):
"""
Returns GCD of multivariate coefficients.

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_content(2*x*y + 6*x + 4*y + 12)
2*y + 6

"""
cont, v = dmp_LC(f, K), u - 1

if dmp_zero_p(f, u):
return cont

for c in f[1:]:
cont = dmp_gcd(cont, c, v, K)

if dmp_one_p(cont, v, K):
break

if K.is_negative(dmp_ground_LC(cont, v, K)):
cont = dmp_neg(cont, v, K)

return cont

[docs]def dmp_primitive(f, u, K):
"""
Returns multivariate content and a primitive polynomial.

Examples
========

>>> R, x, y = ring("x y", ZZ)

>>> R.dmp_primitive(2*x*y + 6*x + 4*y + 12)
(2*y + 6, x + 2)

"""
cont, v = dmp_content(f, u, K), u - 1

if dmp_zero_p(f, u) or dmp_one_p(cont, v, K):
return cont, f
else:
return cont, [dmp_quo(c, cont, v, K) for c in f]