Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
"""Euclidean algorithms, GCDs, LCMs and polynomial remainder sequences. """
from __future__ import print_function, division
from sympy.polys.densebasic import ( dup_strip, dmp_raise, dmp_zero, dmp_one, dmp_ground, dmp_one_p, dmp_zero_p, dmp_zeros, dup_degree, dmp_degree, dmp_degree_in, dup_LC, dmp_LC, dmp_ground_LC, dmp_multi_deflate, dmp_inflate, dup_convert, dmp_convert, dmp_apply_pairs)
from sympy.polys.densearith import ( dup_sub_mul, dup_neg, dmp_neg, dmp_add, dmp_sub, dup_mul, dmp_mul, dmp_pow, dup_div, dmp_div, dup_rem, dup_quo, dmp_quo, dup_prem, dmp_prem, dup_mul_ground, dmp_mul_ground, dmp_mul_term, dup_quo_ground, dmp_quo_ground, dup_max_norm, dmp_max_norm)
from sympy.polys.densetools import ( dup_clear_denoms, dmp_clear_denoms, dup_diff, dmp_diff, dup_eval, dmp_eval, dmp_eval_in, dup_trunc, dmp_ground_trunc, dup_monic, dmp_ground_monic, dup_primitive, dmp_ground_primitive, dup_extract, dmp_ground_extract)
from sympy.polys.galoistools import ( gf_int, gf_crt)
from sympy.polys.polyerrors import ( MultivariatePolynomialError, HeuristicGCDFailed, HomomorphismFailed, NotInvertible, DomainError)
from sympy.polys.polyconfig import query
from sympy.ntheory import nextprime
from sympy.core.compatibility import range
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 ========
>>> from sympy.polys import ring, QQ >>> 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.has_Field: raise DomainError("can't compute half extended GCD over %s" % K)
a, b = [K.one], []
while g: q, r = dup_div(f, g, K) f, g = g, r a, b = b, dup_sub_mul(a, q, b, K)
a = dup_quo_ground(a, dup_LC(f, K), K) f = dup_monic(f, K)
return a, f
def dmp_half_gcdex(f, g, u, K): """ Half extended Euclidean algorithm in `F[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x,y = ring("x,y", ZZ)
""" if not u: return dup_half_gcdex(f, g, K) else: raise MultivariatePolynomialError(f, g)
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 ========
>>> from sympy.polys import ring, QQ >>> 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 = dup_sub_mul(h, s, f, K) t = dup_quo(F, g, K)
return s, t, h
def dmp_gcdex(f, g, u, K): """ Extended Euclidean algorithm in `F[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x,y = ring("x,y", ZZ)
""" if not u: return dup_gcdex(f, g, K) else: raise MultivariatePolynomialError(f, g)
def dup_invert(f, g, K): """ Compute multiplicative inverse of `f` modulo `g` in `F[x]`.
Examples ========
>>> from sympy.polys import ring, QQ >>> 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 dup_rem(s, g, K) else: raise NotInvertible("zero divisor")
def dmp_invert(f, g, u, K): """ Compute multiplicative inverse of `f` modulo `g` in `F[X]`.
Examples ========
>>> from sympy.polys import ring, QQ >>> R, x = ring("x", QQ)
""" if not u: return dup_invert(f, g, K) else: raise MultivariatePolynomialError(f, g)
def dup_euclidean_prs(f, g, K): """ Euclidean polynomial remainder sequence (PRS) in `K[x]`.
Examples ========
>>> from sympy.polys import ring, QQ >>> 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 = dup_rem(f, g, K)
while h: prs.append(h) f, g = g, h h = dup_rem(f, g, K)
return prs
def dmp_euclidean_prs(f, g, u, K): """ Euclidean polynomial remainder sequence (PRS) in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x,y = ring("x,y", ZZ)
""" if not u: return dup_euclidean_prs(f, g, K) else: raise MultivariatePolynomialError(f, g)
def dup_primitive_prs(f, g, K): """ Primitive polynomial remainder sequence (PRS) in `K[x]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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 = dup_primitive(dup_prem(f, g, K), K)
while h: prs.append(h) f, g = g, h _, h = dup_primitive(dup_prem(f, g, K), K)
return prs
def dmp_primitive_prs(f, g, u, K): """ Primitive polynomial remainder sequence (PRS) in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x,y = ring("x,y", ZZ)
""" if not u: return dup_primitive_prs(f, g, K) else: raise MultivariatePolynomialError(f, g)
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 ========
>>> from sympy.polys import ring, ZZ >>> 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 ==========
[1] W.S. Brown, The Subresultant PRS Algorithm. ACM Transaction of Mathematical Software 4 (1978) 237-249
"""
f, g = g, f n, m = m, n
return [], []
return [f], [K.one]
# Conventional first scalar subdeterminant is 1
q = c**(d - 1) c = K.quo((-lc)**d, q) else:
def dup_subresultants(f, g, K): """ Computes subresultant PRS of two polynomials in `K[x]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_subresultants(x**2 + 1, x**2 - 1) [x**2 + 1, x**2 - 1, -2]
""" return dup_inner_subresultants(f, g, K)[0]
def dup_prs_resultant(f, g, K): """ Resultant algorithm in `K[x]` using subresultant PRS.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_prs_resultant(x**2 + 1, x**2 - 1) (4, [x**2 + 1, x**2 - 1, -2])
""" return (K.zero, [])
return (K.zero, R)
def dup_resultant(f, g, K, includePRS=False): """ Computes resultant of two polynomials in `K[x]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_resultant(x**2 + 1, x**2 - 1) 4
""" return dup_prs_resultant(f, g, K)
def dmp_inner_subresultants(f, g, u, K): """ Subresultant PRS algorithm in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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
""" return dup_inner_subresultants(f, g, K)
f, g = g, f n, m = m, n
return [], []
return [f], [dmp_ground(K.one, v)]
dmp_pow(c, d, v, K), v, K)
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:
def dmp_subresultants(f, g, u, K): """ Computes subresultant PRS of two polynomials in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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]
def dmp_prs_resultant(f, g, u, K): """ Resultant algorithm in `K[X]` using subresultant PRS.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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
""" return dup_prs_resultant(f, g, K)
return (dmp_zero(u - 1), [])
return (dmp_zero(u - 1), R)
def dmp_zz_modular_resultant(f, g, p, u, K): """ Compute resultant of `f` and `g` modulo a prime `p`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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(f, u) m = dmp_degree(g, 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 dup_degree(D) <= 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(F, v) == n: G = dmp_eval_in(g, gf_int(a, p), 1, u, K)
if dmp_degree(G, v) == m: break
R = dmp_zz_modular_resultant(F, G, p, v, K) e = dmp_eval(r, a, v, K)
if not v: R = dup_strip([R]) e = dup_strip([e]) else: R = [R] e = [e]
d = K.invert(dup_eval(D, a, K), p) d = dup_mul_ground(D, d, 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)
def dmp_zz_collins_resultant(f, g, u, K): """ Collins's modular resultant algorithm in `Z[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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(f, u) m = dmp_degree(g, 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 K.is_one(P): r = R else: r = dmp_apply_pairs(r, R, _collins_crt, (P, p, K), v, K)
P *= p
return r
def dmp_qq_collins_resultant(f, g, u, K0): """ Collins's modular resultant algorithm in `Q[X]`.
Examples ========
>>> from sympy.polys import ring, QQ >>> R, x,y = ring("x,y", QQ)
>>> f = QQ(1,2)*x + 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(f, u) m = dmp_degree(g, u)
if n < 0 or m < 0: return dmp_zero(u - 1)
K1 = K0.get_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)
def dmp_resultant(f, g, u, K, includePRS=False): """ Computes resultant of two polynomials in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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
"""
return dmp_prs_resultant(f, g, u, K)
if K.is_QQ and query('USE_COLLINS_RESULTANT'): return dmp_qq_collins_resultant(f, g, u, K) else: return dmp_zz_collins_resultant(f, g, u, K)
def dup_discriminant(f, K): """ Computes discriminant of a polynomial in `K[x]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_discriminant(x**2 + 2*x + 3) -8
""" d = dup_degree(f)
if d <= 0: return K.zero else: s = (-1)**((d*(d - 1)) // 2) c = dup_LC(f, K)
r = dup_resultant(f, dup_diff(f, 1, K), K)
return K.quo(r, c*K(s))
def dmp_discriminant(f, u, K): """ Computes discriminant of a polynomial in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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
""" if not u: return dup_discriminant(f, K)
d, v = dmp_degree(f, u), u - 1
if d <= 0: return dmp_zero(v) else: s = (-1)**((d*(d - 1)) // 2) c = dmp_LC(f, K)
r = dmp_resultant(f, dmp_diff(f, 1, u, K), u, K) c = dmp_mul_ground(c, K(s), v, K)
return dmp_quo(r, c, v, K)
def _dup_rr_trivial_gcd(f, g, K): """Handle trivial cases in GCD algorithm over a ring. """ return [], [], [] else: else:
def _dup_ff_trivial_gcd(f, g, K): """Handle trivial cases in GCD algorithm over a field. """ if not (f or g): return [], [], [] elif not f: return dup_monic(g, K), [], [dup_LC(g, K)] elif not g: return dup_monic(f, K), [dup_LC(f, K)], [] else: return None
def _dmp_rr_trivial_gcd(f, g, u, K): """Handle trivial cases in GCD algorithm over a ring. """
return tuple(dmp_zeros(3, u, K)) else: else: else: return None
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 query('USE_SIMPLIFY_GCD'): return _dmp_simplify_gcd(f, g, u, K) else: return None
def _dmp_simplify_gcd(f, g, u, K): """Try to eliminate `x_0` from GCD computation in `K[X]`. """
else: else:
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 ========
>>> from sympy.polys import ring, ZZ >>> 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 = _dup_rr_trivial_gcd(f, g, K)
if result is not None: return result
fc, F = dup_primitive(f, K) gc, G = dup_primitive(g, K)
c = K.gcd(fc, gc)
h = dup_subresultants(F, G, K)[-1] _, h = dup_primitive(h, K)
if K.is_negative(dup_LC(h, K)): c = -c
h = dup_mul_ground(h, c, K)
cff = dup_quo(f, h, K) cfg = dup_quo(g, h, 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 ========
>>> from sympy.polys import ring, QQ >>> 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 = _dup_ff_trivial_gcd(f, g, K)
if result is not None: return result
h = dup_subresultants(f, g, K)[-1] h = dup_monic(h, K)
cff = dup_quo(f, h, K) cfg = dup_quo(g, h, K)
return h, cff, cfg
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 ========
>>> from sympy.polys import ring, ZZ >>> 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
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 ========
>>> from sympy.polys import ring, QQ >>> R, x,y, = ring("x,y", QQ)
>>> f = QQ(1,2)*x**2 + x*y + QQ(1,2)*y**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
HEU_GCD_MAX = 6
def _dup_zz_gcd_interpolate(h, x, K): """Interpolate polynomial GCD from integer GCD. """
def dup_zz_heu_gcd(f, g, K): """ Heuristic polynomial GCD in `Z[x]`.
Given univariate polynomials `f` and `g` in `Z[x]`, returns their GCD and cofactors, i.e. polynomials ``h``, ``cff`` and ``cfg`` such that::
h = gcd(f, g), cff = quo(f, h) and cfg = quo(g, h)
The algorithm is purely heuristic which means it may fail to compute the GCD. This will be signaled by raising an exception. In this case you will need to switch to another GCD method.
The algorithm computes the polynomial GCD by evaluating polynomials f and g at certain points and computing (fast) integer GCD of those evaluations. The polynomial GCD is recovered from the integer image by interpolation. The final step is to verify if the result is the correct GCD. This gives cofactors as a side effect.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_zz_heu_gcd(x**2 - 1, x**2 - 3*x + 2) (x - 1, x + 1, x - 2)
References ==========
1. [Liao95]_
"""
2*min(f_norm // abs(dup_LC(f, K)), g_norm // abs(dup_LC(g, K))) + 2)
h = dup_mul_ground(h, gcd, K) return h, cff, cfg_
cff_, r = dup_div(f, h, K)
if not r: h = dup_mul_ground(h, gcd, K) return h, cff_, cfg
raise HeuristicGCDFailed('no luck')
def _dmp_zz_gcd_interpolate(h, x, v, K): """Interpolate polynomial GCD from integer GCD. """
return dmp_neg(f, v + 1, K) else:
def dmp_zz_heu_gcd(f, g, u, K): """ Heuristic polynomial GCD in `Z[X]`.
Given univariate polynomials `f` and `g` in `Z[X]`, returns their GCD and cofactors, i.e. polynomials ``h``, ``cff`` and ``cfg`` such that::
h = gcd(f, g), cff = quo(f, h) and cfg = quo(g, h)
The algorithm is purely heuristic which means it may fail to compute the GCD. This will be signaled by raising an exception. In this case you will need to switch to another GCD method.
The algorithm computes the polynomial GCD by evaluating polynomials f and g at certain points and computing (fast) integer GCD of those evaluations. The polynomial GCD is recovered from the integer image by interpolation. The evaluation proces reduces f and g variable by variable into a large integer. The final step is to verify if the interpolated polynomial is the correct GCD. This gives cofactors of the input polynomials as a side effect.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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)
References ==========
1. [Liao95]_
"""
2*min(f_norm // abs(dmp_ground_LC(f, u, K)), g_norm // abs(dmp_ground_LC(g, u, K))) + 2)
cff = _dmp_zz_gcd_interpolate(cff, x, v, K)
h, r = dmp_div(f, cff, u, K)
if dmp_zero_p(r, u): cfg_, r = dmp_div(g, h, u, K)
if dmp_zero_p(r, u): h = dmp_mul_ground(h, gcd, u, K) return h, cff, cfg_
cfg = _dmp_zz_gcd_interpolate(cfg, x, v, K)
h, r = dmp_div(g, cfg, u, K)
if dmp_zero_p(r, u): cff_, r = dmp_div(f, h, u, K)
if dmp_zero_p(r, u): h = dmp_mul_ground(h, gcd, u, K) return h, cff_, cfg
x = 73794*x * K.sqrt(K.sqrt(x)) // 27011
raise HeuristicGCDFailed('no luck')
def dup_qq_heu_gcd(f, g, 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 ========
>>> from sympy.polys import ring, QQ >>> R, x = ring("x", QQ)
>>> f = QQ(1,2)*x**2 + QQ(7,4)*x + QQ(3,2) >>> g = QQ(1,2)*x**2 + x
>>> R.dup_qq_heu_gcd(f, g) (x + 2, 1/2*x + 3/4, 1/2*x)
""" result = _dup_ff_trivial_gcd(f, g, K0)
if result is not None: return result
K1 = K0.get_ring()
cf, f = dup_clear_denoms(f, K0, K1) cg, g = dup_clear_denoms(g, K0, K1)
f = dup_convert(f, K0, K1) g = dup_convert(g, K0, K1)
h, cff, cfg = dup_zz_heu_gcd(f, g, K1)
h = dup_convert(h, K1, K0)
c = dup_LC(h, K0) h = dup_monic(h, K0)
cff = dup_convert(cff, K1, K0) cfg = dup_convert(cfg, K1, K0)
cff = dup_mul_ground(cff, K0.quo(c, cf), K0) cfg = dup_mul_ground(cfg, K0.quo(c, cg), K0)
return h, cff, cfg
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 ========
>>> from sympy.polys import ring, QQ >>> R, x,y, = ring("x,y", QQ)
>>> f = QQ(1,4)*x**2 + x*y + y**2 >>> g = QQ(1,2)*x**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.get_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 dup_inner_gcd(f, g, K): """ Computes polynomial GCD and cofactors of `f` and `g` in `K[x]`.
Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``, ``cff = quo(f, h)``, and ``cfg = quo(g, h)``.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_inner_gcd(x**2 - 1, x**2 - 3*x + 2) (x - 1, x + 1, x - 2)
""" try: exact = K.get_exact() except DomainError: return [K.one], f, g
f = dup_convert(f, K, exact) g = dup_convert(g, K, exact)
h, cff, cfg = dup_inner_gcd(f, g, exact)
h = dup_convert(h, exact, K) cff = dup_convert(cff, exact, K) cfg = dup_convert(cfg, exact, K)
return h, cff, cfg if K.is_QQ and query('USE_HEU_GCD'): try: return dup_qq_heu_gcd(f, g, K) except HeuristicGCDFailed: pass
return dup_ff_prs_gcd(f, g, K) else: except HeuristicGCDFailed: pass
return dup_rr_prs_gcd(f, g, K)
def _dmp_inner_gcd(f, g, u, K): """Helper function for `dmp_inner_gcd()`. """ 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 if K.is_QQ and query('USE_HEU_GCD'): try: return dmp_qq_heu_gcd(f, g, u, K) except HeuristicGCDFailed: pass
return dmp_ff_prs_gcd(f, g, u, K) else: except HeuristicGCDFailed: pass
return dmp_rr_prs_gcd(f, g, u, K)
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 ``a = gcd(f, g)``, ``cff = quo(f, h)``, and ``cfg = quo(g, h)``.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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)
"""
dmp_inflate(cff, J, u, K), dmp_inflate(cfg, J, u, K))
def dup_gcd(f, g, K): """ Computes polynomial GCD of `f` and `g` in `K[x]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_gcd(x**2 - 1, x**2 - 3*x + 2) x - 1
"""
def dmp_gcd(f, g, u, K): """ Computes polynomial GCD of `f` and `g` in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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
"""
def dup_rr_lcm(f, g, K): """ Computes polynomial LCM over a ring in `K[x]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_rr_lcm(x**2 - 1, x**2 - 3*x + 2) x**3 - 2*x**2 - x + 2
""" fc, f = dup_primitive(f, K) gc, g = dup_primitive(g, K)
c = K.lcm(fc, gc)
h = dup_quo(dup_mul(f, g, K), dup_gcd(f, g, K), K)
return dup_mul_ground(h, c, K)
def dup_ff_lcm(f, g, K): """ Computes polynomial LCM over a field in `K[x]`.
Examples ========
>>> from sympy.polys import ring, QQ >>> R, x = ring("x", QQ)
>>> f = QQ(1,2)*x**2 + QQ(7,4)*x + QQ(3,2) >>> g = QQ(1,2)*x**2 + x
>>> R.dup_ff_lcm(f, g) x**3 + 7/2*x**2 + 3*x
""" h = dup_quo(dup_mul(f, g, K), dup_gcd(f, g, K), K)
return dup_monic(h, K)
def dup_lcm(f, g, K): """ Computes polynomial LCM of `f` and `g` in `K[x]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_lcm(x**2 - 1, x**2 - 3*x + 2) x**3 - 2*x**2 - x + 2
""" if K.has_Field: return dup_ff_lcm(f, g, K) else: return dup_rr_lcm(f, g, K)
def dmp_rr_lcm(f, g, u, K): """ Computes polynomial LCM over a ring in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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 ========
>>> from sympy.polys import ring, QQ >>> R, x,y, = ring("x,y", QQ)
>>> f = QQ(1,4)*x**2 + x*y + y**2 >>> g = QQ(1,2)*x**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)
def dmp_lcm(f, g, u, K): """ Computes polynomial LCM of `f` and `g` in `K[X]`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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 not u: return dup_lcm(f, g, K)
if K.has_Field: return dmp_ff_lcm(f, g, u, K) else: return dmp_rr_lcm(f, g, u, K)
def dmp_content(f, u, K): """ Returns GCD of multivariate coefficients.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x,y, = ring("x,y", ZZ)
>>> R.dmp_content(2*x*y + 6*x + 4*y + 12) 2*y + 6
"""
return cont
return dmp_neg(cont, v, K) else:
def dmp_primitive(f, u, K): """ Returns multivariate content and a primitive polynomial.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x,y, = ring("x,y", ZZ)
>>> R.dmp_primitive(2*x*y + 6*x + 4*y + 12) (2*y + 6, x + 2)
"""
else:
def dup_cancel(f, g, K, include=True): """ Cancel common factors in a rational function `f/g`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ)
>>> R.dup_cancel(2*x**2 - 2, x**2 - 2*x + 1) (2*x + 2, x - 1)
""" return dmp_cancel(f, g, 0, K, include=include)
def dmp_cancel(f, g, u, K, include=True): """ Cancel common factors in a rational function `f/g`.
Examples ========
>>> from sympy.polys import ring, ZZ >>> 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, K = K, K.get_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.cofactors(cp, cq)
p = dmp_convert(p, u, K, K0) q = dmp_convert(q, u, K, K0)
K = K0
p, q = dmp_neg(p, u, K), dmp_neg(q, u, K)
p = dmp_mul_ground(p, cp, u, K) q = dmp_mul_ground(q, cq, u, K)
return p, q |