from sympy.core.symbol import Dummy from sympy.ntheory import nextprime from sympy.ntheory.modular import crt from sympy.polys.domains import PolynomialRing from sympy.polys.galoistools import ( gf_gcd, gf_from_dict, gf_gcdex, gf_div, gf_lcm) from sympy.polys.polyerrors import ModularGCDFailed from mpmath import sqrt import random def _trivial_gcd(f, g): """ Compute the GCD of two polynomials in trivial cases, i.e. when one or both polynomials are zero. """ ring = f.ring if not (f or g): return ring.zero, ring.zero, ring.zero elif not f: if g.LC < ring.domain.zero: return -g, ring.zero, -ring.one else: return g, ring.zero, ring.one elif not g: if f.LC < ring.domain.zero: return -f, -ring.one, ring.zero else: return f, ring.one, ring.zero return None def _gf_gcd(fp, gp, p): r""" Compute the GCD of two univariate polynomials in `\mathbb{Z}_p[x]`. """ dom = fp.ring.domain while gp: rem = fp deg = gp.degree() lcinv = dom.invert(gp.LC, p) while True: degrem = rem.degree() if degrem < deg: break rem = (rem - gp.mul_monom((degrem - deg,)).mul_ground(lcinv * rem.LC)).trunc_ground(p) fp = gp gp = rem return fp.mul_ground(dom.invert(fp.LC, p)).trunc_ground(p) def _degree_bound_univariate(f, g): r""" Compute an upper bound for the degree of the GCD of two univariate integer polynomials `f` and `g`. The function chooses a suitable prime `p` and computes the GCD of `f` and `g` in `\mathbb{Z}_p[x]`. The choice of `p` guarantees that the degree in `\mathbb{Z}_p[x]` is greater than or equal to the degree in `\mathbb{Z}[x]`. Parameters ========== f : PolyElement univariate integer polynomial g : PolyElement univariate integer polynomial """ gamma = f.ring.domain.gcd(f.LC, g.LC) p = 1 p = nextprime(p) while gamma % p == 0: p = nextprime(p) fp = f.trunc_ground(p) gp = g.trunc_ground(p) hp = _gf_gcd(fp, gp, p) deghp = hp.degree() return deghp def _chinese_remainder_reconstruction_univariate(hp, hq, p, q): r""" Construct a polynomial `h_{pq}` in `\mathbb{Z}_{p q}[x]` such that .. math :: h_{pq} = h_p \; \mathrm{mod} \, p h_{pq} = h_q \; \mathrm{mod} \, q for relatively prime integers `p` and `q` and polynomials `h_p` and `h_q` in `\mathbb{Z}_p[x]` and `\mathbb{Z}_q[x]` respectively. The coefficients of the polynomial `h_{pq}` are computed with the Chinese Remainder Theorem. The symmetric representation in `\mathbb{Z}_p[x]`, `\mathbb{Z}_q[x]` and `\mathbb{Z}_{p q}[x]` is used. It is assumed that `h_p` and `h_q` have the same degree. Parameters ========== hp : PolyElement univariate integer polynomial with coefficients in `\mathbb{Z}_p` hq : PolyElement univariate integer polynomial with coefficients in `\mathbb{Z}_q` p : Integer modulus of `h_p`, relatively prime to `q` q : Integer modulus of `h_q`, relatively prime to `p` Examples ======== >>> from sympy.polys.modulargcd import _chinese_remainder_reconstruction_univariate >>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ) >>> p = 3 >>> q = 5 >>> hp = -x**3 - 1 >>> hq = 2*x**3 - 2*x**2 + x >>> hpq = _chinese_remainder_reconstruction_univariate(hp, hq, p, q) >>> hpq 2*x**3 + 3*x**2 + 6*x + 5 >>> hpq.trunc_ground(p) == hp True >>> hpq.trunc_ground(q) == hq True """ n = hp.degree() x = hp.ring.gens[0] hpq = hp.ring.zero for i in range(n+1): hpq[(i,)] = crt([p, q], [hp.coeff(x**i), hq.coeff(x**i)], symmetric=True)[0] hpq.strip_zero() return hpq def modgcd_univariate(f, g): r""" Computes the GCD of two polynomials in `\mathbb{Z}[x]` using a modular algorithm. The algorithm computes the GCD of two univariate integer polynomials `f` and `g` by computing the GCD in `\mathbb{Z}_p[x]` for suitable primes `p` and then reconstructing the coefficients with the Chinese Remainder Theorem. Trial division is only made for candidates which are very likely the desired GCD. Parameters ========== f : PolyElement univariate integer polynomial g : PolyElement univariate integer polynomial Returns ======= h : PolyElement GCD of the polynomials `f` and `g` cff : PolyElement cofactor of `f`, i.e. `\frac{f}{h}` cfg : PolyElement cofactor of `g`, i.e. `\frac{g}{h}` Examples ======== >>> from sympy.polys.modulargcd import modgcd_univariate >>> from sympy.polys import ring, ZZ >>> R, x = ring("x", ZZ) >>> f = x**5 - 1 >>> g = x - 1 >>> h, cff, cfg = modgcd_univariate(f, g) >>> h, cff, cfg (x - 1, x**4 + x**3 + x**2 + x + 1, 1) >>> cff * h == f True >>> cfg * h == g True >>> f = 6*x**2 - 6 >>> g = 2*x**2 + 4*x + 2 >>> h, cff, cfg = modgcd_univariate(f, g) >>> h, cff, cfg (2*x + 2, 3*x - 3, x + 1) >>> cff * h == f True >>> cfg * h == g True References ========== 1. [Monagan00]_ """ assert f.ring == g.ring and f.ring.domain.is_ZZ result = _trivial_gcd(f, g) if result is not None: return result ring = f.ring cf, f = f.primitive() cg, g = g.primitive() ch = ring.domain.gcd(cf, cg) bound = _degree_bound_univariate(f, g) if bound == 0: return ring(ch), f.mul_ground(cf // ch), g.mul_ground(cg // ch) gamma = ring.domain.gcd(f.LC, g.LC) m = 1 p = 1 while True: p = nextprime(p) while gamma % p == 0: p = nextprime(p) fp = f.trunc_ground(p) gp = g.trunc_ground(p) hp = _gf_gcd(fp, gp, p) deghp = hp.degree() if deghp > bound: continue elif deghp < bound: m = 1 bound = deghp continue hp = hp.mul_ground(gamma).trunc_ground(p) if m == 1: m = p hlastm = hp continue hm = _chinese_remainder_reconstruction_univariate(hp, hlastm, p, m) m *= p if not hm == hlastm: hlastm = hm continue h = hm.quo_ground(hm.content()) fquo, frem = f.div(h) gquo, grem = g.div(h) if not frem and not grem: if h.LC < 0: ch = -ch h = h.mul_ground(ch) cff = fquo.mul_ground(cf // ch) cfg = gquo.mul_ground(cg // ch) return h, cff, cfg def _primitive(f, p): r""" Compute the content and the primitive part of a polynomial in `\mathbb{Z}_p[x_0, \ldots, x_{k-2}, y] \cong \mathbb{Z}_p[y][x_0, \ldots, x_{k-2}]`. Parameters ========== f : PolyElement integer polynomial in `\mathbb{Z}_p[x0, \ldots, x{k-2}, y]` p : Integer modulus of `f` Returns ======= contf : PolyElement integer polynomial in `\mathbb{Z}_p[y]`, content of `f` ppf : PolyElement primitive part of `f`, i.e. `\frac{f}{contf}` Examples ======== >>> from sympy.polys.modulargcd import _primitive >>> from sympy.polys import ring, ZZ >>> R, x, y = ring("x, y", ZZ) >>> p = 3 >>> f = x**2*y**2 + x**2*y - y**2 - y >>> _primitive(f, p) (y**2 + y, x**2 - 1) >>> R, x, y, z = ring("x, y, z", ZZ) >>> f = x*y*z - y**2*z**2 >>> _primitive(f, p) (z, x*y - y**2*z) """ ring = f.ring dom = ring.domain k = ring.ngens coeffs = {} for monom, coeff in f.iterterms(): if monom[:-1] not in coeffs: coeffs[monom[:-1]] = {} coeffs[monom[:-1]][monom[-1]] = coeff cont = [] for coeff in iter(coeffs.values()): cont = gf_gcd(cont, gf_from_dict(coeff, p, dom), p, dom) yring = ring.clone(symbols=ring.symbols[k-1]) contf = yring.from_dense(cont).trunc_ground(p) return contf, f.quo(contf.set_ring(ring)) def _deg(f): r""" Compute the degree of a multivariate polynomial `f \in K[x_0, \ldots, x_{k-2}, y] \cong K[y][x_0, \ldots, x_{k-2}]`. Parameters ========== f : PolyElement polynomial in `K[x_0, \ldots, x_{k-2}, y]` Returns ======= degf : Integer tuple degree of `f` in `x_0, \ldots, x_{k-2}` Examples ======== >>> from sympy.polys.modulargcd import _deg >>> from sympy.polys import ring, ZZ >>> R, x, y = ring("x, y", ZZ) >>> f = x**2*y**2 + x**2*y - 1 >>> _deg(f) (2,) >>> R, x, y, z = ring("x, y, z", ZZ) >>> f = x**2*y**2 + x**2*y - 1 >>> _deg(f) (2, 2) >>> f = x*y*z - y**2*z**2 >>> _deg(f) (1, 1) """ k = f.ring.ngens degf = (0,) * (k-1) for monom in f.itermonoms(): if monom[:-1] > degf: degf = monom[:-1] return degf def _LC(f): r""" Compute the leading coefficient of a multivariate polynomial `f \in K[x_0, \ldots, x_{k-2}, y] \cong K[y][x_0, \ldots, x_{k-2}]`. Parameters ========== f : PolyElement polynomial in `K[x_0, \ldots, x_{k-2}, y]` Returns ======= lcf : PolyElement polynomial in `K[y]`, leading coefficient of `f` Examples ======== >>> from sympy.polys.modulargcd import _LC >>> from sympy.polys import ring, ZZ >>> R, x, y = ring("x, y", ZZ) >>> f = x**2*y**2 + x**2*y - 1 >>> _LC(f) y**2 + y >>> R, x, y, z = ring("x, y, z", ZZ) >>> f = x**2*y**2 + x**2*y - 1 >>> _LC(f) 1 >>> f = x*y*z - y**2*z**2 >>> _LC(f) z """ ring = f.ring k = ring.ngens yring = ring.clone(symbols=ring.symbols[k-1]) y = yring.gens[0] degf = _deg(f) lcf = yring.zero for monom, coeff in f.iterterms(): if monom[:-1] == degf: lcf += coeff*y**monom[-1] return lcf def _swap(f, i): """ Make the variable `x_i` the leading one in a multivariate polynomial `f`. """ ring = f.ring fswap = ring.zero for monom, coeff in f.iterterms(): monomswap = (monom[i],) + monom[:i] + monom[i+1:] fswap[monomswap] = coeff return fswap def _degree_bound_bivariate(f, g): r""" Compute upper degree bounds for the GCD of two bivariate integer polynomials `f` and `g`. The GCD is viewed as a polynomial in `\mathbb{Z}[y][x]` and the function returns an upper bound for its degree and one for the degree of its content. This is done by choosing a suitable prime `p` and computing the GCD of the contents of `f \; \mathrm{mod} \, p` and `g \; \mathrm{mod} \, p`. The choice of `p` guarantees that the degree of the content in `\mathbb{Z}_p[y]` is greater than or equal to the degree in `\mathbb{Z}[y]`. To obtain the degree bound in the variable `x`, the polynomials are evaluated at `y = a` for a suitable `a \in \mathbb{Z}_p` and then their GCD in `\mathbb{Z}_p[x]` is computed. If no such `a` exists, i.e. the degree in `\mathbb{Z}_p[x]` is always smaller than the one in `\mathbb{Z}[y][x]`, then the bound is set to the minimum of the degrees of `f` and `g` in `x`. Parameters ========== f : PolyElement bivariate integer polynomial g : PolyElement bivariate integer polynomial Returns ======= xbound : Integer upper bound for the degree of the GCD of the polynomials `f` and `g` in the variable `x` ycontbound : Integer upper bound for the degree of the content of the GCD of the polynomials `f` and `g` in the variable `y` References ========== 1. [Monagan00]_ """ ring = f.ring gamma1 = ring.domain.gcd(f.LC, g.LC) gamma2 = ring.domain.gcd(_swap(f, 1).LC, _swap(g, 1).LC) badprimes = gamma1 * gamma2 p = 1 p = nextprime(p) while badprimes % p == 0: p = nextprime(p) fp = f.trunc_ground(p) gp = g.trunc_ground(p) contfp, fp = _primitive(fp, p) contgp, gp = _primitive(gp, p) conthp = _gf_gcd(contfp, contgp, p) # polynomial in Z_p[y] ycontbound = conthp.degree() # polynomial in Z_p[y] delta = _gf_gcd(_LC(fp), _LC(gp), p) for a in range(p): if not delta.evaluate(0, a) % p: continue fpa = fp.evaluate(1, a).trunc_ground(p) gpa = gp.evaluate(1, a).trunc_ground(p) hpa = _gf_gcd(fpa, gpa, p) xbound = hpa.degree() return xbound, ycontbound return min(fp.degree(), gp.degree()), ycontbound def _chinese_remainder_reconstruction_multivariate(hp, hq, p, q): r""" Construct a polynomial `h_{pq}` in `\mathbb{Z}_{p q}[x_0, \ldots, x_{k-1}]` such that .. math :: h_{pq} = h_p \; \mathrm{mod} \, p h_{pq} = h_q \; \mathrm{mod} \, q for relatively prime integers `p` and `q` and polynomials `h_p` and `h_q` in `\mathbb{Z}_p[x_0, \ldots, x_{k-1}]` and `\mathbb{Z}_q[x_0, \ldots, x_{k-1}]` respectively. The coefficients of the polynomial `h_{pq}` are computed with the Chinese Remainder Theorem. The symmetric representation in `\mathbb{Z}_p[x_0, \ldots, x_{k-1}]`, `\mathbb{Z}_q[x_0, \ldots, x_{k-1}]` and `\mathbb{Z}_{p q}[x_0, \ldots, x_{k-1}]` is used. Parameters ========== hp : PolyElement multivariate integer polynomial with coefficients in `\mathbb{Z}_p` hq : PolyElement multivariate integer polynomial with coefficients in `\mathbb{Z}_q` p : Integer modulus of `h_p`, relatively prime to `q` q : Integer modulus of `h_q`, relatively prime to `p` Examples ======== >>> from sympy.polys.modulargcd import _chinese_remainder_reconstruction_multivariate >>> from sympy.polys import ring, ZZ >>> R, x, y = ring("x, y", ZZ) >>> p = 3 >>> q = 5 >>> hp = x**3*y - x**2 - 1 >>> hq = -x**3*y - 2*x*y**2 + 2 >>> hpq = _chinese_remainder_reconstruction_multivariate(hp, hq, p, q) >>> hpq 4*x**3*y + 5*x**2 + 3*x*y**2 + 2 >>> hpq.trunc_ground(p) == hp True >>> hpq.trunc_ground(q) == hq True >>> R, x, y, z = ring("x, y, z", ZZ) >>> p = 6 >>> q = 5 >>> hp = 3*x**4 - y**3*z + z >>> hq = -2*x**4 + z >>> hpq = _chinese_remainder_reconstruction_multivariate(hp, hq, p, q) >>> hpq 3*x**4 + 5*y**3*z + z >>> hpq.trunc_ground(p) == hp True >>> hpq.trunc_ground(q) == hq True """ hpmonoms = set(hp.monoms()) hqmonoms = set(hq.monoms()) monoms = hpmonoms.intersection(hqmonoms) hpmonoms.difference_update(monoms) hqmonoms.difference_update(monoms) zero = hp.ring.domain.zero hpq = hp.ring.zero if isinstance(hp.ring.domain, PolynomialRing): crt_ = _chinese_remainder_reconstruction_multivariate else: def crt_(cp, cq, p, q): return crt([p, q], [cp, cq], symmetric=True)[0] for monom in monoms: hpq[monom] = crt_(hp[monom], hq[monom], p, q) for monom in hpmonoms: hpq[monom] = crt_(hp[monom], zero, p, q) for monom in hqmonoms: hpq[monom] = crt_(zero, hq[monom], p, q) return hpq def _interpolate_multivariate(evalpoints, hpeval, ring, i, p, ground=False): r""" Reconstruct a polynomial `h_p` in `\mathbb{Z}_p[x_0, \ldots, x_{k-1}]` from a list of evaluation points in `\mathbb{Z}_p` and a list of polynomials in `\mathbb{Z}_p[x_0, \ldots, x_{i-1}, x_{i+1}, \ldots, x_{k-1}]`, which are the images of `h_p` evaluated in the variable `x_i`. It is also possible to reconstruct a parameter of the ground domain, i.e. if `h_p` is a polynomial over `\mathbb{Z}_p[x_0, \ldots, x_{k-1}]`. In this case, one has to set ``ground=True``. Parameters ========== evalpoints : list of Integer objects list of evaluation points in `\mathbb{Z}_p` hpeval : list of PolyElement objects list of polynomials in (resp. over) `\mathbb{Z}_p[x_0, \ldots, x_{i-1}, x_{i+1}, \ldots, x_{k-1}]`, images of `h_p` evaluated in the variable `x_i` ring : PolyRing `h_p` will be an element of this ring i : Integer index of the variable which has to be reconstructed p : Integer prime number, modulus of `h_p` ground : Boolean indicates whether `x_i` is in the ground domain, default is ``False`` Returns ======= hp : PolyElement interpolated polynomial in (resp. over) `\mathbb{Z}_p[x_0, \ldots, x_{k-1}]` """ hp = ring.zero if ground: domain = ring.domain.domain y = ring.domain.gens[i] else: domain = ring.domain y = ring.gens[i] for a, hpa in zip(evalpoints, hpeval): numer = ring.one denom = domain.one for b in evalpoints: if b == a: continue numer *= y - b denom *= a - b denom = domain.invert(denom, p) coeff = numer.mul_ground(denom) hp += hpa.set_ring(ring) * coeff return hp.trunc_ground(p) def modgcd_bivariate(f, g): r""" Computes the GCD of two polynomials in `\mathbb{Z}[x, y]` using a modular algorithm. The algorithm computes the GCD of two bivariate integer polynomials `f` and `g` by calculating the GCD in `\mathbb{Z}_p[x, y]` for suitable primes `p` and then reconstructing the coefficients with the Chinese Remainder Theorem. To compute the bivariate GCD over `\mathbb{Z}_p`, the polynomials `f \; \mathrm{mod} \, p` and `g \; \mathrm{mod} \, p` are evaluated at `y = a` for certain `a \in \mathbb{Z}_p` and then their univariate GCD in `\mathbb{Z}_p[x]` is computed. Interpolating those yields the bivariate GCD in `\mathbb{Z}_p[x, y]`. To verify the result in `\mathbb{Z}[x, y]`, trial division is done, but only for candidates which are very likely the desired GCD. Parameters ========== f : PolyElement bivariate integer polynomial g : PolyElement bivariate integer polynomial Returns ======= h : PolyElement GCD of the polynomials `f` and `g` cff : PolyElement cofactor of `f`, i.e. `\frac{f}{h}` cfg : PolyElement cofactor of `g`, i.e. `\frac{g}{h}` Examples ======== >>> from sympy.polys.modulargcd import modgcd_bivariate >>> from sympy.polys import ring, ZZ >>> R, x, y = ring("x, y", ZZ) >>> f = x**2 - y**2 >>> g = x**2 + 2*x*y + y**2 >>> h, cff, cfg = modgcd_bivariate(f, g) >>> h, cff, cfg (x + y, x - y, x + y) >>> cff * h == f True >>> cfg * h == g True >>> f = x**2*y - x**2 - 4*y + 4 >>> g = x + 2 >>> h, cff, cfg = modgcd_bivariate(f, g) >>> h, cff, cfg (x + 2, x*y - x - 2*y + 2, 1) >>> cff * h == f True >>> cfg * h == g True References ========== 1. [Monagan00]_ """ assert f.ring == g.ring and f.ring.domain.is_ZZ result = _trivial_gcd(f, g) if result is not None: return result ring = f.ring cf, f = f.primitive() cg, g = g.primitive() ch = ring.domain.gcd(cf, cg) xbound, ycontbound = _degree_bound_bivariate(f, g) if xbound == ycontbound == 0: return ring(ch), f.mul_ground(cf // ch), g.mul_ground(cg // ch) fswap = _swap(f, 1) gswap = _swap(g, 1) degyf = fswap.degree() degyg = gswap.degree() ybound, xcontbound = _degree_bound_bivariate(fswap, gswap) if ybound == xcontbound == 0: return ring(ch), f.mul_ground(cf // ch), g.mul_ground(cg // ch) # TODO: to improve performance, choose the main variable here gamma1 = ring.domain.gcd(f.LC, g.LC) gamma2 = ring.domain.gcd(fswap.LC, gswap.LC) badprimes = gamma1 * gamma2 m = 1 p = 1 while True: p = nextprime(p) while badprimes % p == 0: p = nextprime(p) fp = f.trunc_ground(p) gp = g.trunc_ground(p) contfp, fp = _primitive(fp, p) contgp, gp = _primitive(gp, p) conthp = _gf_gcd(contfp, contgp, p) # monic polynomial in Z_p[y] degconthp = conthp.degree() if degconthp > ycontbound: continue elif degconthp < ycontbound: m = 1 ycontbound = degconthp continue # polynomial in Z_p[y] delta = _gf_gcd(_LC(fp), _LC(gp), p) degcontfp = contfp.degree() degcontgp = contgp.degree() degdelta = delta.degree() N = min(degyf - degcontfp, degyg - degcontgp, ybound - ycontbound + degdelta) + 1 if p < N: continue n = 0 evalpoints = [] hpeval = [] unlucky = False for a in range(p): deltaa = delta.evaluate(0, a) if not deltaa % p: continue fpa = fp.evaluate(1, a).trunc_ground(p) gpa = gp.evaluate(1, a).trunc_ground(p) hpa = _gf_gcd(fpa, gpa, p) # monic polynomial in Z_p[x] deghpa = hpa.degree() if deghpa > xbound: continue elif deghpa < xbound: m = 1 xbound = deghpa unlucky = True break hpa = hpa.mul_ground(deltaa).trunc_ground(p) evalpoints.append(a) hpeval.append(hpa) n += 1 if n == N: break if unlucky: continue if n < N: continue hp = _interpolate_multivariate(evalpoints, hpeval, ring, 1, p) hp = _primitive(hp, p)[1] hp = hp * conthp.set_ring(ring) degyhp = hp.degree(1) if degyhp > ybound: continue if degyhp < ybound: m = 1 ybound = degyhp continue hp = hp.mul_ground(gamma1).trunc_ground(p) if m == 1: m = p hlastm = hp continue hm = _chinese_remainder_reconstruction_multivariate(hp, hlastm, p, m) m *= p if not hm == hlastm: hlastm = hm continue h = hm.quo_ground(hm.content()) fquo, frem = f.div(h) gquo, grem = g.div(h) if not frem and not grem: if h.LC < 0: ch = -ch h = h.mul_ground(ch) cff = fquo.mul_ground(cf // ch) cfg = gquo.mul_ground(cg // ch) return h, cff, cfg def _modgcd_multivariate_p(f, g, p, degbound, contbound): r""" Compute the GCD of two polynomials in `\mathbb{Z}_p[x_0, \ldots, x_{k-1}]`. The algorithm reduces the problem step by step by evaluating the polynomials `f` and `g` at `x_{k-1} = a` for suitable `a \in \mathbb{Z}_p` and then calls itself recursively to compute the GCD in `\mathbb{Z}_p[x_0, \ldots, x_{k-2}]`. If these recursive calls are successful for enough evaluation points, the GCD in `k` variables is interpolated, otherwise the algorithm returns ``None``. Every time a GCD or a content is computed, their degrees are compared with the bounds. If a degree greater then the bound is encountered, then the current call returns ``None`` and a new evaluation point has to be chosen. If at some point the degree is smaller, the correspondent bound is updated and the algorithm fails. Parameters ========== f : PolyElement multivariate integer polynomial with coefficients in `\mathbb{Z}_p` g : PolyElement multivariate integer polynomial with coefficients in `\mathbb{Z}_p` p : Integer prime number, modulus of `f` and `g` degbound : list of Integer objects ``degbound[i]`` is an upper bound for the degree of the GCD of `f` and `g` in the variable `x_i` contbound : list of Integer objects ``contbound[i]`` is an upper bound for the degree of the content of the GCD in `\mathbb{Z}_p[x_i][x_0, \ldots, x_{i-1}]`, ``contbound[0]`` is not used can therefore be chosen arbitrarily. Returns ======= h : PolyElement GCD of the polynomials `f` and `g` or ``None`` References ========== 1. [Monagan00]_ 2. [Brown71]_ """ ring = f.ring k = ring.ngens if k == 1: h = _gf_gcd(f, g, p).trunc_ground(p) degh = h.degree() if degh > degbound[0]: return None if degh < degbound[0]: degbound[0] = degh raise ModularGCDFailed return h degyf = f.degree(k-1) degyg = g.degree(k-1) contf, f = _primitive(f, p) contg, g = _primitive(g, p) conth = _gf_gcd(contf, contg, p) # polynomial in Z_p[y] degcontf = contf.degree() degcontg = contg.degree() degconth = conth.degree() if degconth > contbound[k-1]: return None if degconth < contbound[k-1]: contbound[k-1] = degconth raise ModularGCDFailed lcf = _LC(f) lcg = _LC(g) delta = _gf_gcd(lcf, lcg, p) # polynomial in Z_p[y] evaltest = delta for i in range(k-1): evaltest *= _gf_gcd(_LC(_swap(f, i)), _LC(_swap(g, i)), p) degdelta = delta.degree() N = min(degyf - degcontf, degyg - degcontg, degbound[k-1] - contbound[k-1] + degdelta) + 1 if p < N: return None n = 0 d = 0 evalpoints = [] heval = [] points = list(range(p)) while points: a = random.sample(points, 1)[0] points.remove(a) if not evaltest.evaluate(0, a) % p: continue deltaa = delta.evaluate(0, a) % p fa = f.evaluate(k-1, a).trunc_ground(p) ga = g.evaluate(k-1, a).trunc_ground(p) # polynomials in Z_p[x_0, ..., x_{k-2}] ha = _modgcd_multivariate_p(fa, ga, p, degbound, contbound) if ha is None: d += 1 if d > n: return None continue if ha.is_ground: h = conth.set_ring(ring).trunc_ground(p) return h ha = ha.mul_ground(deltaa).trunc_ground(p) evalpoints.append(a) heval.append(ha) n += 1 if n == N: h = _interpolate_multivariate(evalpoints, heval, ring, k-1, p) h = _primitive(h, p)[1] * conth.set_ring(ring) degyh = h.degree(k-1) if degyh > degbound[k-1]: return None if degyh < degbound[k-1]: degbound[k-1] = degyh raise ModularGCDFailed return h return None def modgcd_multivariate(f, g): r""" Compute the GCD of two polynomials in `\mathbb{Z}[x_0, \ldots, x_{k-1}]` using a modular algorithm. The algorithm computes the GCD of two multivariate integer polynomials `f` and `g` by calculating the GCD in `\mathbb{Z}_p[x_0, \ldots, x_{k-1}]` for suitable primes `p` and then reconstructing the coefficients with the Chinese Remainder Theorem. To compute the multivariate GCD over `\mathbb{Z}_p` the recursive subroutine :func:`_modgcd_multivariate_p` is used. To verify the result in `\mathbb{Z}[x_0, \ldots, x_{k-1}]`, trial division is done, but only for candidates which are very likely the desired GCD. Parameters ========== f : PolyElement multivariate integer polynomial g : PolyElement multivariate integer polynomial Returns ======= h : PolyElement GCD of the polynomials `f` and `g` cff : PolyElement cofactor of `f`, i.e. `\frac{f}{h}` cfg : PolyElement cofactor of `g`, i.e. `\frac{g}{h}` Examples ======== >>> from sympy.polys.modulargcd import modgcd_multivariate >>> from sympy.polys import ring, ZZ >>> R, x, y = ring("x, y", ZZ) >>> f = x**2 - y**2 >>> g = x**2 + 2*x*y + y**2 >>> h, cff, cfg = modgcd_multivariate(f, g) >>> h, cff, cfg (x + y, x - y, x + y) >>> cff * h == f True >>> cfg * h == g True >>> R, x, y, z = ring("x, y, z", ZZ) >>> f = x*z**2 - y*z**2 >>> g = x**2*z + z >>> h, cff, cfg = modgcd_multivariate(f, g) >>> h, cff, cfg (z, x*z - y*z, x**2 + 1) >>> cff * h == f True >>> cfg * h == g True References ========== 1. [Monagan00]_ 2. [Brown71]_ See also ======== _modgcd_multivariate_p """ assert f.ring == g.ring and f.ring.domain.is_ZZ result = _trivial_gcd(f, g) if result is not None: return result ring = f.ring k = ring.ngens # divide out integer content cf, f = f.primitive() cg, g = g.primitive() ch = ring.domain.gcd(cf, cg) gamma = ring.domain.gcd(f.LC, g.LC) badprimes = ring.domain.one for i in range(k): badprimes *= ring.domain.gcd(_swap(f, i).LC, _swap(g, i).LC) degbound = [min(fdeg, gdeg) for fdeg, gdeg in zip(f.degrees(), g.degrees())] contbound = list(degbound) m = 1 p = 1 while True: p = nextprime(p) while badprimes % p == 0: p = nextprime(p) fp = f.trunc_ground(p) gp = g.trunc_ground(p) try: # monic GCD of fp, gp in Z_p[x_0, ..., x_{k-2}, y] hp = _modgcd_multivariate_p(fp, gp, p, degbound, contbound) except ModularGCDFailed: m = 1 continue if hp is None: continue hp = hp.mul_ground(gamma).trunc_ground(p) if m == 1: m = p hlastm = hp continue hm = _chinese_remainder_reconstruction_multivariate(hp, hlastm, p, m) m *= p if not hm == hlastm: hlastm = hm continue h = hm.primitive()[1] fquo, frem = f.div(h) gquo, grem = g.div(h) if not frem and not grem: if h.LC < 0: ch = -ch h = h.mul_ground(ch) cff = fquo.mul_ground(cf // ch) cfg = gquo.mul_ground(cg // ch) return h, cff, cfg def _gf_div(f, g, p): r""" Compute `\frac f g` modulo `p` for two univariate polynomials over `\mathbb Z_p`. """ ring = f.ring densequo, denserem = gf_div(f.to_dense(), g.to_dense(), p, ring.domain) return ring.from_dense(densequo), ring.from_dense(denserem) def _rational_function_reconstruction(c, p, m): r""" Reconstruct a rational function `\frac a b` in `\mathbb Z_p(t)` from .. math:: c = \frac a b \; \mathrm{mod} \, m, where `c` and `m` are polynomials in `\mathbb Z_p[t]` and `m` has positive degree. The algorithm is based on the Euclidean Algorithm. In general, `m` is not irreducible, so it is possible that `b` is not invertible modulo `m`. In that case ``None`` is returned. Parameters ========== c : PolyElement univariate polynomial in `\mathbb Z[t]` p : Integer prime number m : PolyElement modulus, not necessarily irreducible Returns ======= frac : FracElement either `\frac a b` in `\mathbb Z(t)` or ``None`` References ========== 1. [Hoeij04]_ """ ring = c.ring domain = ring.domain M = m.degree() N = M // 2 D = M - N - 1 r0, s0 = m, ring.zero r1, s1 = c, ring.one while r1.degree() > N: quo = _gf_div(r0, r1, p)[0] r0, r1 = r1, (r0 - quo*r1).trunc_ground(p) s0, s1 = s1, (s0 - quo*s1).trunc_ground(p) a, b = r1, s1 if b.degree() > D or _gf_gcd(b, m, p) != 1: return None lc = b.LC if lc != 1: lcinv = domain.invert(lc, p) a = a.mul_ground(lcinv).trunc_ground(p) b = b.mul_ground(lcinv).trunc_ground(p) field = ring.to_field() return field(a) / field(b) def _rational_reconstruction_func_coeffs(hm, p, m, ring, k): r""" Reconstruct every coefficient `c_h` of a polynomial `h` in `\mathbb Z_p(t_k)[t_1, \ldots, t_{k-1}][x, z]` from the corresponding coefficient `c_{h_m}` of a polynomial `h_m` in `\mathbb Z_p[t_1, \ldots, t_k][x, z] \cong \mathbb Z_p[t_k][t_1, \ldots, t_{k-1}][x, z]` such that .. math:: c_{h_m} = c_h \; \mathrm{mod} \, m, where `m \in \mathbb Z_p[t]`. The reconstruction is based on the Euclidean Algorithm. In general, `m` is not irreducible, so it is possible that this fails for some coefficient. In that case ``None`` is returned. Parameters ========== hm : PolyElement polynomial in `\mathbb Z[t_1, \ldots, t_k][x, z]` p : Integer prime number, modulus of `\mathbb Z_p` m : PolyElement modulus, polynomial in `\mathbb Z[t]`, not necessarily irreducible ring : PolyRing `\mathbb Z(t_k)[t_1, \ldots, t_{k-1}][x, z]`, `h` will be an element of this ring k : Integer index of the parameter `t_k` which will be reconstructed Returns ======= h : PolyElement reconstructed polynomial in `\mathbb Z(t_k)[t_1, \ldots, t_{k-1}][x, z]` or ``None`` See also ======== _rational_function_reconstruction """ h = ring.zero for monom, coeff in hm.iterterms(): if k == 0: coeffh = _rational_function_reconstruction(coeff, p, m) if not coeffh: return None else: coeffh = ring.domain.zero for mon, c in coeff.drop_to_ground(k).iterterms(): ch = _rational_function_reconstruction(c, p, m) if not ch: return None coeffh[mon] = ch h[monom] = coeffh return h def _gf_gcdex(f, g, p): r""" Extended Euclidean Algorithm for two univariate polynomials over `\mathbb Z_p`. Returns polynomials `s, t` and `h`, such that `h` is the GCD of `f` and `g` and `sf + tg = h \; \mathrm{mod} \, p`. """ ring = f.ring s, t, h = gf_gcdex(f.to_dense(), g.to_dense(), p, ring.domain) return ring.from_dense(s), ring.from_dense(t), ring.from_dense(h) def _trunc(f, minpoly, p): r""" Compute the reduced representation of a polynomial `f` in `\mathbb Z_p[z] / (\check m_{\alpha}(z))[x]` Parameters ========== f : PolyElement polynomial in `\mathbb Z[x, z]` minpoly : PolyElement polynomial `\check m_{\alpha} \in \mathbb Z[z]`, not necessarily irreducible p : Integer prime number, modulus of `\mathbb Z_p` Returns ======= ftrunc : PolyElement polynomial in `\mathbb Z[x, z]`, reduced modulo `\check m_{\alpha}(z)` and `p` """ ring = f.ring minpoly = minpoly.set_ring(ring) p_ = ring.ground_new(p) return f.trunc_ground(p).rem([minpoly, p_]).trunc_ground(p) def _euclidean_algorithm(f, g, minpoly, p): r""" Compute the monic GCD of two univariate polynomials in `\mathbb{Z}_p[z]/(\check m_{\alpha}(z))[x]` with the Euclidean Algorithm. In general, `\check m_{\alpha}(z)` is not irreducible, so it is possible that some leading coefficient is not invertible modulo `\check m_{\alpha}(z)`. In that case ``None`` is returned. Parameters ========== f, g : PolyElement polynomials in `\mathbb Z[x, z]` minpoly : PolyElement polynomial in `\mathbb Z[z]`, not necessarily irreducible p : Integer prime number, modulus of `\mathbb Z_p` Returns ======= h : PolyElement GCD of `f` and `g` in `\mathbb Z[z, x]` or ``None``, coefficients are in `\left[ -\frac{p-1} 2, \frac{p-1} 2 \right]` """ ring = f.ring f = _trunc(f, minpoly, p) g = _trunc(g, minpoly, p) while g: rem = f deg = g.degree(0) # degree in x lcinv, _, gcd = _gf_gcdex(ring.dmp_LC(g), minpoly, p) if not gcd == 1: return None while True: degrem = rem.degree(0) # degree in x if degrem < deg: break quo = (lcinv * ring.dmp_LC(rem)).set_ring(ring) rem = _trunc(rem - g.mul_monom((degrem - deg, 0))*quo, minpoly, p) f = g g = rem lcfinv = _gf_gcdex(ring.dmp_LC(f), minpoly, p)[0].set_ring(ring) return _trunc(f * lcfinv, minpoly, p) def _trial_division(f, h, minpoly, p=None): r""" Check if `h` divides `f` in `\mathbb K[t_1, \ldots, t_k][z]/(m_{\alpha}(z))`, where `\mathbb K` is either `\mathbb Q` or `\mathbb Z_p`. This algorithm is based on pseudo division and does not use any fractions. By default `\mathbb K` is `\mathbb Q`, if a prime number `p` is given, `\mathbb Z_p` is chosen instead. Parameters ========== f, h : PolyElement polynomials in `\mathbb Z[t_1, \ldots, t_k][x, z]` minpoly : PolyElement polynomial `m_{\alpha}(z)` in `\mathbb Z[t_1, \ldots, t_k][z]` p : Integer or None if `p` is given, `\mathbb K` is set to `\mathbb Z_p` instead of `\mathbb Q`, default is ``None`` Returns ======= rem : PolyElement remainder of `\frac f h` References ========== .. [1] [Hoeij02]_ """ ring = f.ring zxring = ring.clone(symbols=(ring.symbols[1], ring.symbols[0])) minpoly = minpoly.set_ring(ring) rem = f degrem = rem.degree() degh = h.degree() degm = minpoly.degree(1) lch = _LC(h).set_ring(ring) lcm = minpoly.LC while rem and degrem >= degh: # polynomial in Z[t_1, ..., t_k][z] lcrem = _LC(rem).set_ring(ring) rem = rem*lch - h.mul_monom((degrem - degh, 0))*lcrem if p: rem = rem.trunc_ground(p) degrem = rem.degree(1) while rem and degrem >= degm: # polynomial in Z[t_1, ..., t_k][x] lcrem = _LC(rem.set_ring(zxring)).set_ring(ring) rem = rem.mul_ground(lcm) - minpoly.mul_monom((0, degrem - degm))*lcrem if p: rem = rem.trunc_ground(p) degrem = rem.degree(1) degrem = rem.degree() return rem def _evaluate_ground(f, i, a): r""" Evaluate a polynomial `f` at `a` in the `i`-th variable of the ground domain. """ ring = f.ring.clone(domain=f.ring.domain.ring.drop(i)) fa = ring.zero for monom, coeff in f.iterterms(): fa[monom] = coeff.evaluate(i, a) return fa def _func_field_modgcd_p(f, g, minpoly, p): r""" Compute the GCD of two polynomials `f` and `g` in `\mathbb Z_p(t_1, \ldots, t_k)[z]/(\check m_\alpha(z))[x]`. The algorithm reduces the problem step by step by evaluating the polynomials `f` and `g` at `t_k = a` for suitable `a \in \mathbb Z_p` and then calls itself recursively to compute the GCD in `\mathbb Z_p(t_1, \ldots, t_{k-1})[z]/(\check m_\alpha(z))[x]`. If these recursive calls are successful, the GCD over `k` variables is interpolated, otherwise the algorithm returns ``None``. After interpolation, Rational Function Reconstruction is used to obtain the correct coefficients. If this fails, a new evaluation point has to be chosen, otherwise the desired polynomial is obtained by clearing denominators. The result is verified with a fraction free trial division. Parameters ========== f, g : PolyElement polynomials in `\mathbb Z[t_1, \ldots, t_k][x, z]` minpoly : PolyElement polynomial in `\mathbb Z[t_1, \ldots, t_k][z]`, not necessarily irreducible p : Integer prime number, modulus of `\mathbb Z_p` Returns ======= h : PolyElement primitive associate in `\mathbb Z[t_1, \ldots, t_k][x, z]` of the GCD of the polynomials `f` and `g` or ``None``, coefficients are in `\left[ -\frac{p-1} 2, \frac{p-1} 2 \right]` References ========== 1. [Hoeij04]_ """ ring = f.ring domain = ring.domain # Z[t_1, ..., t_k] if isinstance(domain, PolynomialRing): k = domain.ngens else: return _euclidean_algorithm(f, g, minpoly, p) if k == 1: qdomain = domain.ring.to_field() else: qdomain = domain.ring.drop_to_ground(k - 1) qdomain = qdomain.clone(domain=qdomain.domain.ring.to_field()) qring = ring.clone(domain=qdomain) # = Z(t_k)[t_1, ..., t_{k-1}][x, z] n = 1 d = 1 # polynomial in Z_p[t_1, ..., t_k][z] gamma = ring.dmp_LC(f) * ring.dmp_LC(g) # polynomial in Z_p[t_1, ..., t_k] delta = minpoly.LC evalpoints = [] heval = [] LMlist = [] points = list(range(p)) while points: a = random.sample(points, 1)[0] points.remove(a) if k == 1: test = delta.evaluate(k-1, a) % p == 0 else: test = delta.evaluate(k-1, a).trunc_ground(p) == 0 if test: continue gammaa = _evaluate_ground(gamma, k-1, a) minpolya = _evaluate_ground(minpoly, k-1, a) if gammaa.rem([minpolya, gammaa.ring(p)]) == 0: continue fa = _evaluate_ground(f, k-1, a) ga = _evaluate_ground(g, k-1, a) # polynomial in Z_p[x, t_1, ..., t_{k-1}, z]/(minpoly) ha = _func_field_modgcd_p(fa, ga, minpolya, p) if ha is None: d += 1 if d > n: return None continue if ha == 1: return ha LM = [ha.degree()] + [0]*(k-1) if k > 1: for monom, coeff in ha.iterterms(): if monom[0] == LM[0] and coeff.LM > tuple(LM[1:]): LM[1:] = coeff.LM evalpoints_a = [a] heval_a = [ha] if k == 1: m = qring.domain.get_ring().one else: m = qring.domain.domain.get_ring().one t = m.ring.gens[0] for b, hb, LMhb in zip(evalpoints, heval, LMlist): if LMhb == LM: evalpoints_a.append(b) heval_a.append(hb) m *= (t - b) m = m.trunc_ground(p) evalpoints.append(a) heval.append(ha) LMlist.append(LM) n += 1 # polynomial in Z_p[t_1, ..., t_k][x, z] h = _interpolate_multivariate(evalpoints_a, heval_a, ring, k-1, p, ground=True) # polynomial in Z_p(t_k)[t_1, ..., t_{k-1}][x, z] h = _rational_reconstruction_func_coeffs(h, p, m, qring, k-1) if h is None: continue if k == 1: dom = qring.domain.field den = dom.ring.one for coeff in h.itercoeffs(): den = dom.ring.from_dense(gf_lcm(den.to_dense(), coeff.denom.to_dense(), p, dom.domain)) else: dom = qring.domain.domain.field den = dom.ring.one for coeff in h.itercoeffs(): for c in coeff.itercoeffs(): den = dom.ring.from_dense(gf_lcm(den.to_dense(), c.denom.to_dense(), p, dom.domain)) den = qring.domain_new(den.trunc_ground(p)) h = ring(h.mul_ground(den).as_expr()).trunc_ground(p) if not _trial_division(f, h, minpoly, p) and not _trial_division(g, h, minpoly, p): return h return None def _integer_rational_reconstruction(c, m, domain): r""" Reconstruct a rational number `\frac a b` from .. math:: c = \frac a b \; \mathrm{mod} \, m, where `c` and `m` are integers. The algorithm is based on the Euclidean Algorithm. In general, `m` is not a prime number, so it is possible that `b` is not invertible modulo `m`. In that case ``None`` is returned. Parameters ========== c : Integer `c = \frac a b \; \mathrm{mod} \, m` m : Integer modulus, not necessarily prime domain : IntegerRing `a, b, c` are elements of ``domain`` Returns ======= frac : Rational either `\frac a b` in `\mathbb Q` or ``None`` References ========== 1. [Wang81]_ """ if c < 0: c += m r0, s0 = m, domain.zero r1, s1 = c, domain.one bound = sqrt(m / 2) # still correct if replaced by ZZ.sqrt(m // 2) ? while r1 >= bound: quo = r0 // r1 r0, r1 = r1, r0 - quo*r1 s0, s1 = s1, s0 - quo*s1 if abs(s1) >= bound: return None if s1 < 0: a, b = -r1, -s1 elif s1 > 0: a, b = r1, s1 else: return None field = domain.get_field() return field(a) / field(b) def _rational_reconstruction_int_coeffs(hm, m, ring): r""" Reconstruct every rational coefficient `c_h` of a polynomial `h` in `\mathbb Q[t_1, \ldots, t_k][x, z]` from the corresponding integer coefficient `c_{h_m}` of a polynomial `h_m` in `\mathbb Z[t_1, \ldots, t_k][x, z]` such that .. math:: c_{h_m} = c_h \; \mathrm{mod} \, m, where `m \in \mathbb Z`. The reconstruction is based on the Euclidean Algorithm. In general, `m` is not a prime number, so it is possible that this fails for some coefficient. In that case ``None`` is returned. Parameters ========== hm : PolyElement polynomial in `\mathbb Z[t_1, \ldots, t_k][x, z]` m : Integer modulus, not necessarily prime ring : PolyRing `\mathbb Q[t_1, \ldots, t_k][x, z]`, `h` will be an element of this ring Returns ======= h : PolyElement reconstructed polynomial in `\mathbb Q[t_1, \ldots, t_k][x, z]` or ``None`` See also ======== _integer_rational_reconstruction """ h = ring.zero if isinstance(ring.domain, PolynomialRing): reconstruction = _rational_reconstruction_int_coeffs domain = ring.domain.ring else: reconstruction = _integer_rational_reconstruction domain = hm.ring.domain for monom, coeff in hm.iterterms(): coeffh = reconstruction(coeff, m, domain) if not coeffh: return None h[monom] = coeffh return h def _func_field_modgcd_m(f, g, minpoly): r""" Compute the GCD of two polynomials in `\mathbb Q(t_1, \ldots, t_k)[z]/(m_{\alpha}(z))[x]` using a modular algorithm. The algorithm computes the GCD of two polynomials `f` and `g` by calculating the GCD in `\mathbb Z_p(t_1, \ldots, t_k)[z] / (\check m_{\alpha}(z))[x]` for suitable primes `p` and the primitive associate `\check m_{\alpha}(z)` of `m_{\alpha}(z)`. Then the coefficients are reconstructed with the Chinese Remainder Theorem and Rational Reconstruction. To compute the GCD over `\mathbb Z_p(t_1, \ldots, t_k)[z] / (\check m_{\alpha})[x]`, the recursive subroutine ``_func_field_modgcd_p`` is used. To verify the result in `\mathbb Q(t_1, \ldots, t_k)[z] / (m_{\alpha}(z))[x]`, a fraction free trial division is used. Parameters ========== f, g : PolyElement polynomials in `\mathbb Z[t_1, \ldots, t_k][x, z]` minpoly : PolyElement irreducible polynomial in `\mathbb Z[t_1, \ldots, t_k][z]` Returns ======= h : PolyElement the primitive associate in `\mathbb Z[t_1, \ldots, t_k][x, z]` of the GCD of `f` and `g` Examples ======== >>> from sympy.polys.modulargcd import _func_field_modgcd_m >>> from sympy.polys import ring, ZZ >>> R, x, z = ring('x, z', ZZ) >>> minpoly = (z**2 - 2).drop(0) >>> f = x**2 + 2*x*z + 2 >>> g = x + z >>> _func_field_modgcd_m(f, g, minpoly) x + z >>> D, t = ring('t', ZZ) >>> R, x, z = ring('x, z', D) >>> minpoly = (z**2-3).drop(0) >>> f = x**2 + (t + 1)*x*z + 3*t >>> g = x*z + 3*t >>> _func_field_modgcd_m(f, g, minpoly) x + t*z References ========== 1. [Hoeij04]_ See also ======== _func_field_modgcd_p """ ring = f.ring domain = ring.domain if isinstance(domain, PolynomialRing): k = domain.ngens QQdomain = domain.ring.clone(domain=domain.domain.get_field()) QQring = ring.clone(domain=QQdomain) else: k = 0 QQring = ring.clone(domain=ring.domain.get_field()) cf, f = f.primitive() cg, g = g.primitive() # polynomial in Z[t_1, ..., t_k][z] gamma = ring.dmp_LC(f) * ring.dmp_LC(g) # polynomial in Z[t_1, ..., t_k] delta = minpoly.LC p = 1 primes = [] hplist = [] LMlist = [] while True: p = nextprime(p) if gamma.trunc_ground(p) == 0: continue if k == 0: test = (delta % p == 0) else: test = (delta.trunc_ground(p) == 0) if test: continue fp = f.trunc_ground(p) gp = g.trunc_ground(p) minpolyp = minpoly.trunc_ground(p) hp = _func_field_modgcd_p(fp, gp, minpolyp, p) if hp is None: continue if hp == 1: return ring.one LM = [hp.degree()] + [0]*k if k > 0: for monom, coeff in hp.iterterms(): if monom[0] == LM[0] and coeff.LM > tuple(LM[1:]): LM[1:] = coeff.LM hm = hp m = p for q, hq, LMhq in zip(primes, hplist, LMlist): if LMhq == LM: hm = _chinese_remainder_reconstruction_multivariate(hq, hm, q, m) m *= q primes.append(p) hplist.append(hp) LMlist.append(LM) hm = _rational_reconstruction_int_coeffs(hm, m, QQring) if hm is None: continue if k == 0: h = hm.clear_denoms()[1] else: den = domain.domain.one for coeff in hm.itercoeffs(): den = domain.domain.lcm(den, coeff.clear_denoms()[0]) h = hm.mul_ground(den) # convert back to Z[t_1, ..., t_k][x, z] from Q[t_1, ..., t_k][x, z] h = h.set_ring(ring) h = h.primitive()[1] if not (_trial_division(f.mul_ground(cf), h, minpoly) or _trial_division(g.mul_ground(cg), h, minpoly)): return h def _to_ZZ_poly(f, ring): r""" Compute an associate of a polynomial `f \in \mathbb Q(\alpha)[x_0, \ldots, x_{n-1}]` in `\mathbb Z[x_1, \ldots, x_{n-1}][z] / (\check m_{\alpha}(z))[x_0]`, where `\check m_{\alpha}(z) \in \mathbb Z[z]` is the primitive associate of the minimal polynomial `m_{\alpha}(z)` of `\alpha` over `\mathbb Q`. Parameters ========== f : PolyElement polynomial in `\mathbb Q(\alpha)[x_0, \ldots, x_{n-1}]` ring : PolyRing `\mathbb Z[x_1, \ldots, x_{n-1}][x_0, z]` Returns ======= f_ : PolyElement associate of `f` in `\mathbb Z[x_1, \ldots, x_{n-1}][x_0, z]` """ f_ = ring.zero if isinstance(ring.domain, PolynomialRing): domain = ring.domain.domain else: domain = ring.domain den = domain.one for coeff in f.itercoeffs(): for c in coeff.rep: if c: den = domain.lcm(den, c.denominator) for monom, coeff in f.iterterms(): coeff = coeff.rep m = ring.domain.one if isinstance(ring.domain, PolynomialRing): m = m.mul_monom(monom[1:]) n = len(coeff) for i in range(n): if coeff[i]: c = domain(coeff[i] * den) * m if (monom[0], n-i-1) not in f_: f_[(monom[0], n-i-1)] = c else: f_[(monom[0], n-i-1)] += c return f_ def _to_ANP_poly(f, ring): r""" Convert a polynomial `f \in \mathbb Z[x_1, \ldots, x_{n-1}][z]/(\check m_{\alpha}(z))[x_0]` to a polynomial in `\mathbb Q(\alpha)[x_0, \ldots, x_{n-1}]`, where `\check m_{\alpha}(z) \in \mathbb Z[z]` is the primitive associate of the minimal polynomial `m_{\alpha}(z)` of `\alpha` over `\mathbb Q`. Parameters ========== f : PolyElement polynomial in `\mathbb Z[x_1, \ldots, x_{n-1}][x_0, z]` ring : PolyRing `\mathbb Q(\alpha)[x_0, \ldots, x_{n-1}]` Returns ======= f_ : PolyElement polynomial in `\mathbb Q(\alpha)[x_0, \ldots, x_{n-1}]` """ domain = ring.domain f_ = ring.zero if isinstance(f.ring.domain, PolynomialRing): for monom, coeff in f.iterterms(): for mon, coef in coeff.iterterms(): m = (monom[0],) + mon c = domain([domain.domain(coef)] + [0]*monom[1]) if m not in f_: f_[m] = c else: f_[m] += c else: for monom, coeff in f.iterterms(): m = (monom[0],) c = domain([domain.domain(coeff)] + [0]*monom[1]) if m not in f_: f_[m] = c else: f_[m] += c return f_ def _minpoly_from_dense(minpoly, ring): r""" Change representation of the minimal polynomial from ``DMP`` to ``PolyElement`` for a given ring. """ minpoly_ = ring.zero for monom, coeff in minpoly.terms(): minpoly_[monom] = ring.domain(coeff) return minpoly_ def _primitive_in_x0(f): r""" Compute the content in `x_0` and the primitive part of a polynomial `f` in `\mathbb Q(\alpha)[x_0, x_1, \ldots, x_{n-1}] \cong \mathbb Q(\alpha)[x_1, \ldots, x_{n-1}][x_0]`. """ fring = f.ring ring = fring.drop_to_ground(*range(1, fring.ngens)) dom = ring.domain.ring f_ = ring(f.as_expr()) cont = dom.zero for coeff in f_.itercoeffs(): cont = func_field_modgcd(cont, coeff)[0] if cont == dom.one: return cont, f return cont, f.quo(cont.set_ring(fring)) # TODO: add support for algebraic function fields def func_field_modgcd(f, g): r""" Compute the GCD of two polynomials `f` and `g` in `\mathbb Q(\alpha)[x_0, \ldots, x_{n-1}]` using a modular algorithm. The algorithm first computes the primitive associate `\check m_{\alpha}(z)` of the minimal polynomial `m_{\alpha}` in `\mathbb{Z}[z]` and the primitive associates of `f` and `g` in `\mathbb{Z}[x_1, \ldots, x_{n-1}][z]/(\check m_{\alpha})[x_0]`. Then it computes the GCD in `\mathbb Q(x_1, \ldots, x_{n-1})[z]/(m_{\alpha}(z))[x_0]`. This is done by calculating the GCD in `\mathbb{Z}_p(x_1, \ldots, x_{n-1})[z]/(\check m_{\alpha}(z))[x_0]` for suitable primes `p` and then reconstructing the coefficients with the Chinese Remainder Theorem and Rational Reconstuction. The GCD over `\mathbb{Z}_p(x_1, \ldots, x_{n-1})[z]/(\check m_{\alpha}(z))[x_0]` is computed with a recursive subroutine, which evaluates the polynomials at `x_{n-1} = a` for suitable evaluation points `a \in \mathbb Z_p` and then calls itself recursively until the ground domain does no longer contain any parameters. For `\mathbb{Z}_p[z]/(\check m_{\alpha}(z))[x_0]` the Euclidean Algorithm is used. The results of those recursive calls are then interpolated and Rational Function Reconstruction is used to obtain the correct coefficients. The results, both in `\mathbb Q(x_1, \ldots, x_{n-1})[z]/(m_{\alpha}(z))[x_0]` and `\mathbb{Z}_p(x_1, \ldots, x_{n-1})[z]/(\check m_{\alpha}(z))[x_0]`, are verified by a fraction free trial division. Apart from the above GCD computation some GCDs in `\mathbb Q(\alpha)[x_1, \ldots, x_{n-1}]` have to be calculated, because treating the polynomials as univariate ones can result in a spurious content of the GCD. For this ``func_field_modgcd`` is called recursively. Parameters ========== f, g : PolyElement polynomials in `\mathbb Q(\alpha)[x_0, \ldots, x_{n-1}]` Returns ======= h : PolyElement monic GCD of the polynomials `f` and `g` cff : PolyElement cofactor of `f`, i.e. `\frac f h` cfg : PolyElement cofactor of `g`, i.e. `\frac g h` Examples ======== >>> from sympy.polys.modulargcd import func_field_modgcd >>> from sympy.polys import AlgebraicField, QQ, ring >>> from sympy import sqrt >>> A = AlgebraicField(QQ, sqrt(2)) >>> R, x = ring('x', A) >>> f = x**2 - 2 >>> g = x + sqrt(2) >>> h, cff, cfg = func_field_modgcd(f, g) >>> h == x + sqrt(2) True >>> cff * h == f True >>> cfg * h == g True >>> R, x, y = ring('x, y', A) >>> f = x**2 + 2*sqrt(2)*x*y + 2*y**2 >>> g = x + sqrt(2)*y >>> h, cff, cfg = func_field_modgcd(f, g) >>> h == x + sqrt(2)*y True >>> cff * h == f True >>> cfg * h == g True >>> f = x + sqrt(2)*y >>> g = x + y >>> h, cff, cfg = func_field_modgcd(f, g) >>> h == R.one True >>> cff * h == f True >>> cfg * h == g True References ========== 1. [Hoeij04]_ """ ring = f.ring domain = ring.domain n = ring.ngens assert ring == g.ring and domain.is_Algebraic result = _trivial_gcd(f, g) if result is not None: return result z = Dummy('z') ZZring = ring.clone(symbols=ring.symbols + (z,), domain=domain.domain.get_ring()) if n == 1: f_ = _to_ZZ_poly(f, ZZring) g_ = _to_ZZ_poly(g, ZZring) minpoly = ZZring.drop(0).from_dense(domain.mod.rep) h = _func_field_modgcd_m(f_, g_, minpoly) h = _to_ANP_poly(h, ring) else: # 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] ZZring_ = ZZring.drop_to_ground(*range(1, n)) f_ = _to_ZZ_poly(f, ZZring_) g_ = _to_ZZ_poly(g, ZZring_) minpoly = _minpoly_from_dense(domain.mod, ZZring_.drop(0)) h = _func_field_modgcd_m(f_, g_, minpoly) h = _to_ANP_poly(h, ring) contx0h_, h = _primitive_in_x0(h) h *= contx0h.set_ring(ring) f *= contx0f.set_ring(ring) g *= contx0g.set_ring(ring) h = h.quo_ground(h.LC) return h, f.quo(h), g.quo(h)