Skip to content

Commit

Permalink
polys: removed _GCD.cofactors(), gcd helpers return only gcd
Browse files Browse the repository at this point in the history
  • Loading branch information
skirpichev committed Apr 5, 2021
1 parent e428c3f commit 072d6e5
Show file tree
Hide file tree
Showing 3 changed files with 59 additions and 96 deletions.
7 changes: 4 additions & 3 deletions diofant/domains/ring.py
Expand Up @@ -51,9 +51,10 @@ def half_gcdex(self, a, b):

def cofactors(self, a, b):
"""Returns GCD and cofactors of ``a`` and ``b``."""
gcd = self.gcd(a, b)
cfa = self.quo(a, gcd)
cfb = self.quo(b, gcd)
gcd, cfa, cfb = self.gcd(a, b), self.zero, self.zero
if gcd:
cfa = self.quo(a, gcd)
cfb = self.quo(b, gcd)
return gcd, cfa, cfb

def lcm(self, a, b):
Expand Down
121 changes: 44 additions & 77 deletions diofant/polys/euclidtools.py
Expand Up @@ -14,7 +14,21 @@ class _GCD:

def gcd(self, f, g):
"""Returns GCD of ``f`` and ``g``."""
return self.cofactors(f, g)[0]
if f.is_zero and g.is_zero:
return self.zero
elif f.is_zero:
return self._gcd_zero(g)
elif g.is_zero:
return self._gcd_zero(f)
elif f.is_term:
return self._gcd_term(f, g)
elif g.is_term:
return self._gcd_term(g, f)

J, (f, g) = self._deflate(f, g)
h = self._gcd(f, g)

return self._inflate(h, J)

def lcm(self, f, g):
"""Returns LCM of ``f`` and ``g``."""
Expand All @@ -32,29 +46,6 @@ def lcm(self, f, g):
else:
return h.monic()

def cofactors(self, f, g):
"""Returns GCD and cofactors of ``f`` and ``g``."""
if f.is_zero and g.is_zero:
zero = self.zero
return zero, zero, zero
elif f.is_zero:
h, cff, cfg = self._gcd_zero(g)
return h, cff, cfg
elif g.is_zero:
h, cfg, cff = self._gcd_zero(f)
return h, cff, cfg
elif f.is_term:
h, cff, cfg = self._gcd_term(f, g)
return h, cff, cfg
elif g.is_term:
h, cfg, cff = self._gcd_term(g, f)
return h, cff, cfg

J, (f, g) = self._deflate(f, g)
h, cff, cfg = self._gcd(f, g)

return self._inflate(h, J), self._inflate(cff, J), self._inflate(cfg, J)

def _deflate(self, *polys):
J = [0]*self.ngens

Expand Down Expand Up @@ -95,21 +86,15 @@ def _inflate(self, f, J):
return poly

def _gcd_zero(self, f):
one, zero = self.one, self.zero
if self.domain.is_Field:
return f.monic(), zero, self.ground_new(f.LC)
return f.monic()
else:
if not self.is_normal(f):
return -f, zero, -one
else:
return f, zero, one
return f if self.is_normal(f) else -f

def _gcd_term(self, f, g):
domain = self.domain
ground_gcd = domain.gcd
ground_quo = domain.quo
mf, cf = f.LT
_mgcd, _cgcd = mf, cf
_mgcd, _cgcd = f.LT
if domain.is_Field:
for mg, cg in g.items():
_mgcd = _mgcd.gcd(mg)
Expand All @@ -118,11 +103,7 @@ def _gcd_term(self, f, g):
for mg, cg in g.items():
_mgcd = _mgcd.gcd(mg)
_cgcd = ground_gcd(_cgcd, cg)
h = self.term_new(_mgcd, _cgcd)
cff = self.term_new(mf/_mgcd, ground_quo(cf, _cgcd))
cfg = self.from_dict({mg/_mgcd: ground_quo(cg, _cgcd)
for mg, cg in g.items()})
return h, cff, cfg
return self.term_new(_mgcd, _cgcd)

def _gcd(self, f, g):
domain = self.domain
Expand All @@ -137,13 +118,13 @@ def _gcd(self, f, g):
try:
exact = domain.get_exact()
except DomainError:
return self.one, f, g
return self.one

f, g = map(operator.methodcaller('set_domain', exact), (f, g))
ring = self.clone(domain=exact)
h = ring._gcd(f, g)

return tuple(map(operator.methodcaller('set_domain', domain),
ring.cofactors(f, g)))
return h.set_domain(domain)
elif domain.is_Field:
return self._ff_prs_gcd(f, g)
else:
Expand All @@ -167,20 +148,15 @@ def _gcd_ZZ(self, f, g):
def _gcd_QQ(self, f, g):
domain = self.domain

cf, f = f.clear_denoms(convert=True)
cg, g = g.clear_denoms(convert=True)
_, f = f.clear_denoms(convert=True)
_, g = g.clear_denoms(convert=True)

ring = self.clone(domain=domain.ring)

h, cff, cfg = map(operator.methodcaller('set_ring', self),
ring._gcd_ZZ(f, g))

c, h = h.LC, h.monic()

cff *= domain.quo(c, cf)
cfg *= domain.quo(c, cg)
h = ring._gcd_ZZ(f, g)
h = h.set_ring(self)

return h, cff, cfg
return h.monic()

def _gcd_AA(self, f, g):
from .modulargcd import func_field_modgcd
Expand All @@ -196,10 +172,9 @@ def _zz_heu_gcd(self, f, g):
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::
their GCD, i.e. polynomial ``h``::
h = gcd(f, g), cff = quo(f, h) and cfg = quo(g, h)
h = gcd(f, g)
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
Expand All @@ -210,8 +185,7 @@ def _zz_heu_gcd(self, f, g):
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.
interpolated polynomial is the correct GCD.
References
==========
Expand All @@ -236,7 +210,7 @@ def _zz_heu_gcd(self, f, g):
2*min(f_norm // abs(f.LC),
g_norm // abs(g.LC)) + 4)

cofactors = domain.cofactors if ring.is_univariate else ring.drop(0)._zz_heu_gcd
cofactors = domain.cofactors if ring.is_univariate else ring.drop(0).cofactors

for i in range(query('HEU_GCD_MAX')):
ff = f.eval(x0, x)
Expand All @@ -247,36 +221,33 @@ def _zz_heu_gcd(self, f, g):
h = ring._gcd_interpolate(h, x)
h = h.primitive()[1]

cff_, r = divmod(f, h)
_, r = divmod(f, h)

if not r:
cfg_, r = divmod(g, h)
_, r = divmod(g, h)

if not r:
h *= gcd
return h, cff_, cfg_
return h*gcd

cff = ring._gcd_interpolate(cff, x)

h, r = divmod(f, cff)

if not r:
cfg_, r = divmod(g, h)
_, r = divmod(g, h)

if not r:
h *= gcd
return h, cff, cfg_
return h*gcd

cfg = ring._gcd_interpolate(cfg, x)

h, r = divmod(g, cfg)

if not r:
cff_, r = divmod(f, h)
_, r = divmod(f, h)

if not r:
h *= gcd
return h, cff_, cfg
return h*gcd

x = 73794*x * domain.sqrt(domain.sqrt(x)) // 27011

Expand Down Expand Up @@ -310,8 +281,8 @@ def _rr_prs_gcd(self, f, g):
if self.is_multivariate:
ring, f, g = map(operator.methodcaller('eject', *self.gens[1:]),
(ring, f, g))
return tuple(map(operator.methodcaller('inject'),
ring._rr_prs_gcd(f, g)))
h = ring._rr_prs_gcd(f, g)
return h.inject()

domain = ring.domain

Expand All @@ -320,11 +291,9 @@ def _rr_prs_gcd(self, f, g):

h = ff.subresultants(fg)[-1]
_, h = h.primitive()

c = domain.gcd(fc, gc)
h *= c

return h, f // h, g // h
return h*c

def _ff_prs_gcd(self, f, g):
"""Computes polynomial GCD using subresultants over a field."""
Expand All @@ -340,19 +309,17 @@ def _ff_prs_gcd(self, f, g):
F, G = map(operator.methodcaller('inject'), (F, G))

h = F.subresultants(G)[-1]
c, _, _ = ring.domain._ff_prs_gcd(fc, gc)
c = ring.domain._ff_prs_gcd(fc, gc)
h = h.eject(*self.gens[1:])
_, h = h.primitive()
h = h.inject()
h *= c
h = h.monic()

return h, f // h, g // h
return h.monic()

h = f.subresultants(g)[-1]
h = h.monic()

return h, f // h, g // h
return h.monic()

def _primitive_prs(self, f, g):
"""
Expand Down
27 changes: 11 additions & 16 deletions diofant/polys/modulargcd.py
Expand Up @@ -384,12 +384,12 @@ def modgcd(f, g):
>>> _, x, y = ring('x y', ZZ)
>>> modgcd((x - y)*(x + y), (x + y)**2)
(x + y, x - y, x + y)
x + y
>>> _, x, y, z = ring('x y z', ZZ)
>>> modgcd((x - y)*z**2, (x**2 + 1)*z)
(z, x*z - y*z, x**2 + 1)
z
References
==========
Expand Down Expand Up @@ -452,13 +452,10 @@ def modgcd(f, g):
continue

h = hm.primitive()[1]
fquo, frem = divmod(f, h)
gquo, grem = divmod(g, h)
_, frem = divmod(f, h)
_, grem = divmod(g, h)
if not frem and not grem:
h *= ch
cff = fquo*(cf // ch)
cfg = gquo*(cg // ch)
return h, cff, cfg
return h*ch


def _rational_function_reconstruction(c, p, m):
Expand Down Expand Up @@ -1303,7 +1300,7 @@ def _primitive_in_x0(f):
cont = dom.zero

for coeff in f_.values():
cont = func_field_modgcd(cont, coeff)[0]
cont = func_field_modgcd(cont, coeff)
if cont == dom.one:
return cont, f

Expand Down Expand Up @@ -1367,15 +1364,15 @@ def func_field_modgcd(f, g):
>>> _, x = ring('x', A)
>>> func_field_modgcd(x**2 - 2, x + sqrt(2))
(x + sqrt(2), x - sqrt(2), 1)
x + sqrt(2)
>>> _, x, y = ring('x y', A)
>>> func_field_modgcd((x + sqrt(2)*y)**2, x + sqrt(2)*y)
(x + sqrt(2)*y, x + sqrt(2)*y, 1)
x + sqrt(2)*y
>>> func_field_modgcd(x + sqrt(2)*y, x + y)
(1, x + sqrt(2)*y, x + y)
1
References
==========
Expand Down Expand Up @@ -1408,7 +1405,7 @@ def func_field_modgcd(f, g):
# contx0f in Q(a)[x_1, ..., x_{n-1}], f in Q(a)[x_0, ..., x_{n-1}]
contx0f, f = _primitive_in_x0(f)
contx0g, g = _primitive_in_x0(g)
contx0h = func_field_modgcd(contx0f, contx0g)[0]
contx0h = func_field_modgcd(contx0f, contx0g)

ZZring_ = ZZring.eject(*range(1, n))

Expand All @@ -1424,6 +1421,4 @@ def func_field_modgcd(f, g):
f *= contx0f.set_ring(ring)
g *= contx0g.set_ring(ring)

h = h.quo_ground(h.LC)

return h, f//h, g//h
return h.quo_ground(h.LC)

0 comments on commit 072d6e5

Please sign in to comment.