# Source code for diofant.polys.euclidtools

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

from ..core import cacheit
from ..ntheory import nextprime
from .densearith import (dmp_add, dmp_div, dmp_max_norm, dmp_mul,
dmp_mul_ground, dmp_mul_term, dmp_neg, dmp_pow,
dmp_prem, dmp_quo, dmp_quo_ground, dmp_rem, dmp_sub,
dmp_sub_mul, dup_mul)
from .densebasic import (dmp_apply_pairs, dmp_convert, dmp_degree_in,
dmp_ground, dmp_ground_LC, dmp_inflate, dmp_LC,
dmp_multi_deflate, dmp_one, dmp_one_p, dmp_raise,
dmp_strip, dmp_zero, dmp_zero_p, dmp_zeros)
from .densetools import (dmp_clear_denoms, dmp_diff_in, dmp_eval_in,
dmp_ground_monic, dmp_ground_primitive,
dmp_ground_trunc, dup_trunc)
from .galoistools import gf_crt, gf_int
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

>>> R.dup_half_gcdex(f, 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

>>> R.dup_gcdex(f, 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 dup_euclidean_prs(f, g, K):
"""
Euclidean polynomial remainder sequence (PRS) in K[x].

Examples
========

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

>>> f = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
>>> g = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21

>>> prs = R.dup_euclidean_prs(f, g)

>>> prs[0]
x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
>>> prs[1]
3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
>>> prs[2]
-5/9*x**4 + 1/9*x**2 - 1/3
>>> prs[3]
-117/25*x**2 - 9*x + 441/25
>>> prs[4]
233150/19773*x - 102500/6591
>>> prs[5]
-1288744821/543589225

"""
prs = [f, g]
h = dmp_rem(f, g, 0, K)

while h:
prs.append(h)
f, g = g, h
h = dmp_rem(f, g, 0, K)

return prs

[docs]def dup_primitive_prs(f, g, K):
"""
Primitive polynomial remainder sequence (PRS) in K[x].

Examples
========

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

>>> f = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
>>> g = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21

>>> prs = R.dup_primitive_prs(f, g)

>>> prs[0]
x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
>>> prs[1]
3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
>>> prs[2]
5*x**4 - x**2 + 3
>>> prs[3]
13*x**2 + 25*x - 49
>>> prs[4]
4663*x - 6150
>>> prs[5]
1

"""
prs = [f, g]
_, h = dmp_ground_primitive(dmp_prem(f, g, 0, K), 0, K)

while h:
prs.append(h)
f, g = g, h
_, h = dmp_ground_primitive(dmp_prem(f, g, 0, K), 0, K)

return prs

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.dup_inner_subresultants(x**2 + 1, x**2 - 1)
([x**2 + 1, x**2 - 1, -2], [1, 1, 4])

References
==========

* [Brown78]_
"""
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)

>>> R.dup_prs_resultant(x**2 + 1, x**2 - 1)
(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

def dup_resultant(f, g, K, includePRS=False):
"""
Computes resultant of two polynomials in K[x].

Examples
========

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

>>> R.dup_resultant(x**2 + 1, x**2 - 1)
4
"""
if includePRS:
return dup_prs_resultant(f, g, K)
return dup_prs_resultant(f, g, K)[0]

[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

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

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

[docs]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_prs_resultant(f, g)

>>> 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 gf_int(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, gf_int(a, p), 1, u, K)

if dmp_degree_in(F, 0, v) == n:
G = dmp_eval_in(g, gf_int(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 = dup_trunc(D, p, K)

return r

def _collins_crt(r, R, P, p, K):
"""Wrapper of CRT for Collins's resultant algorithm. """
return gf_int(gf_crt([r, R], [P, p], K), P*p)

[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

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

"""
if not u:
return dup_resultant(f, g, K, includePRS=includePRS)

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]

[docs]def dmp_discriminant(f, u, K):
"""
Computes discriminant of a polynomial in K[X].

Examples
========

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

>>> R.dmp_discriminant(x**2*y + x*z + t)
-4*y*t + z**2
"""
d, v = dmp_degree_in(f, 0, u), u - 1

if d <= 0:
return dmp_zero(v) if u else K.zero
else:
s = (-1)**((d*(d - 1)) // 2)
c = dmp_LC(f, K)

r = dmp_resultant(f, dmp_diff_in(f, 1, 0, u, K), u, K)

if u:
c = dmp_mul_ground(c, K(s), v, K)
return dmp_quo(r, c, v, K)
else:
return K.quo(r, c*K(s))

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_nonnegative(dmp_ground_LC(g, u, K)):
return g, dmp_zero(u), dmp_one(u, K)
else:
return dmp_neg(g, u, K), dmp_zero(u), dmp_ground(-K.one, u)
elif zero_g:
if K.is_nonnegative(dmp_ground_LC(f, u, K)):
return f, dmp_one(u, K), dmp_zero(u)
else:
return dmp_neg(f, u, K), dmp_ground(-K.one, u), 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)

========

diofant.polys.heuristicgcd.heugcd

References
==========

* [Liao95]_
"""
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))
h, cff, cfg = heugcd(f, g)
return tuple(map(ring.to_dense, f.cofactors(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))
h, cff, cfg = modgcd(f, g)
return tuple(map(ring.to_dense, f.cofactors(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))
h, cff, cfg = func_field_modgcd(f, g)
return tuple(map(ring.to_dense, f.cofactors(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_multi_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]

def dmp_rr_lcm(f, g, u, K):
"""
Computes polynomial LCM over a ring in K[X].

Examples
========

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

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

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

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

>>> R.dmp_rr_lcm(f, g)
x**3 + 2*x**2*y + x*y**2

"""
fc, f = dmp_ground_primitive(f, u, K)
gc, g = dmp_ground_primitive(g, u, K)

c = K.lcm(fc, gc)

h = dmp_quo(dmp_mul(f, g, u, K),
dmp_gcd(f, g, u, K), u, K)

return dmp_mul_ground(h, c, u, K)

def dmp_ff_lcm(f, g, u, K):
"""
Computes polynomial LCM over a field in K[X].

Examples
========

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

>>> f = (x**2 + 7*x/2 + 3)/2
>>> g = x**2/2 + x

>>> R.dmp_ff_lcm(f, g)
x**3 + 7/2*x**2 + 3*x

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

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

>>> R.dmp_ff_lcm(f, g)
x**3 + 4*x**2*y + 4*x*y**2
"""
h = dmp_quo(dmp_mul(f, g, u, K),
dmp_gcd(f, g, u, K), u, K)

return dmp_ground_monic(h, u, K)

[docs]def dmp_lcm(f, g, u, K):
"""
Computes polynomial LCM 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_lcm(f, g)
x**3 + 2*x**2*y + x*y**2
"""
if K.is_Field:
return dmp_ff_lcm(f, g, u, K)
else:
return dmp_rr_lcm(f, g, u, K)

[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)):
return dmp_neg(cont, v, K)
else:
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 ]

[docs]def dmp_cancel(f, g, u, K, include=True):
"""
Cancel common factors in a rational function f/g.

Examples
========

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

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

"""
K0 = None

if K.is_Field and K.has_assoc_Ring:
K0, K = K, K.ring

cq, f = dmp_clear_denoms(f, u, K0, K, convert=True)
cp, g = dmp_clear_denoms(g, u, K0, K, convert=True)
else:
cp, cq = K.one, K.one

_, p, q = dmp_inner_gcd(f, g, u, K)

if K0 is not None:
_, cp, cq = K.cofactors(cp, cq)

p = dmp_convert(p, u, K, K0)
q = dmp_convert(q, u, K, K0)

K = K0

p_neg = K.is_negative(dmp_ground_LC(p, u, K))
q_neg = K.is_negative(dmp_ground_LC(q, u, K))

if p_neg and q_neg:
p, q = dmp_neg(p, u, K), dmp_neg(q, u, K)
elif p_neg:
cp, p = -cp, dmp_neg(p, u, K)
elif q_neg:
cp, q = -cp, dmp_neg(q, u, K)

if not include:
return cp, cq, p, q

p = dmp_mul_ground(p, cp, u, K)
q = dmp_mul_ground(q, cq, u, K)

return p, q