/

# Source code for sympy.polys.distributedpolys

"""Sparse distributed multivariate polynomials. """

from sympy.polys.monomialtools import (
monomial_mul, monomial_div, monomial_lcm, lex,
)

from sympy.polys.polyerrors import (
ExactQuotientFailed,
)

[docs]def sdp_LC(f, K):
"""Returns the leading coeffcient of f. """
if not f:
return K.zero
else:
return f[0][1]

[docs]def sdp_LM(f, u):
"""Returns the leading monomial of f. """
if not f:
return (0,) * (u + 1)
else:
return f[0][0]

[docs]def sdp_LT(f, u, K):
"""Returns the leading term of f. """
if f:
return f[0]
else:
return (0,) * (u + 1), K.zero

[docs]def sdp_del_LT(f):
"""Removes the leading from f. """
return f[1:]

def sdp_coeffs(f):
"""Returns a list of monomials in f. """
return [ coeff for _, coeff in f ]

[docs]def sdp_monoms(f):
"""Returns a list of monomials in f. """
return [ monom for monom, _ in f ]

[docs]def sdp_sort(f, O):
"""Sort terms in f using the given monomial order O. """
return sorted(f, key=lambda term: O(term[0]), reverse=True)

[docs]def sdp_strip(f):
"""Remove terms with zero coefficients from f in K[X]. """
return [ (monom, coeff) for monom, coeff in f if coeff ]

[docs]def sdp_normal(f, K):
"""Normalize distributed polynomial in the given domain. """
return [ (monom, K.convert(coeff)) for monom, coeff in f if coeff ]

[docs]def sdp_from_dict(f, O):
"""Make a distributed polynomial from a dictionary. """
return sdp_sort(list(f.items()), O)

[docs]def sdp_to_dict(f):
"""Make a dictionary from a distributed polynomial. """
return dict(f)

[docs]def sdp_indep_p(f, j, u):
"""Returns True if a polynomial is independent of x_j. """
if j < 0 or j > u:
raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))
else:
return all(not monom[j] for monom in sdp_monoms(h))

[docs]def sdp_one_p(f, u, K):
"""Returns True if f is a multivariate one in K[X]. """
return f == sdp_one(u, K)

[docs]def sdp_one(u, K):
"""Returns a multivariate one in K[X]. """
return (((0,) * (u + 1), K.one),)

[docs]def sdp_term_p(f):
"""Returns True if f has a single term or is zero. """
return len(f) <= 1

[docs]def sdp_abs(f, u, O, K):
"""Make all coefficients positive in K[X]. """
return [ (monom, K.abs(coeff)) for monom, coeff in f ]

[docs]def sdp_neg(f, u, O, K):
"""Negate a polynomial in K[X]. """
return [ (monom, -coeff) for monom, coeff in f ]

[docs]def sdp_add_term(f, term, u, O, K):
"""Add a single term using bisection method. """
M, c = term

if not c:
return f
if not f:
return [(M, c)]

monoms = sdp_monoms(f)

if O(M) > O(monoms[ 0]):
return [(M, c)] + f
if O(M) < O(monoms[-1]):
return f + [(M, c)]

lo, hi = 0, len(monoms) - 1

while lo <= hi:
i = (lo + hi) // 2

if O(M) == O(monoms[i]):
coeff = f[i][1] + c

if not coeff:
return f[:i] + f[i + 1:]
else:
return f[:i] + [(M, coeff)] + f[i + 1:]
else:
if O(M) > O(monoms[i]):
hi = i - 1
else:
lo = i + 1
else:
return f[:i] + [(M, c)] + f[i + 1:]

[docs]def sdp_sub_term(f, term, u, O, K):
"""Sub a single term using bisection method. """
M, c = term

if not c:
return f
if not f:
return [(M, -c)]

monoms = sdp_monoms(f)

if O(M) > O(monoms[ 0]):
return [(M, -c)] + f
if O(M) < O(monoms[-1]):
return f + [(M, -c)]

lo, hi = 0, len(monoms) - 1

while lo <= hi:
i = (lo + hi) // 2

if O(M) == O(monoms[i]):
coeff = f[i][1] - c

if not coeff:
return f[:i] + f[i + 1:]
else:
return f[:i] + [(M, coeff)] + f[i + 1:]
else:
if O(M) > O(monoms[i]):
hi = i - 1
else:
lo = i + 1
else:
return f[:i] + [(M, -c)] + f[i + 1:]

[docs]def sdp_mul_term(f, term, u, O, K):
"""Multiply a distributed polynomial by a term. """
M, c = term

if not f or not c:
return []
else:
if K.is_one(c):
return [ (monomial_mul(f_M, M), f_c) for f_M, f_c in f ]
else:
return [ (monomial_mul(f_M, M), f_c * c) for f_M, f_c in f ]

[docs]def sdp_add(f, g, u, O, K):
"""Add distributed polynomials in K[X]. """
h = dict(f)

for monom, c in g:
if monom in h:
coeff = h[monom] + c

if not coeff:
del h[monom]
else:
h[monom] = coeff
else:
h[monom] = c

return sdp_from_dict(h, O)

[docs]def sdp_sub(f, g, u, O, K):
"""Subtract distributed polynomials in K[X]. """
h = dict(f)

for monom, c in g:
if monom in h:
coeff = h[monom] - c

if not coeff:
del h[monom]
else:
h[monom] = coeff
else:
h[monom] = -c

return sdp_from_dict(h, O)

[docs]def sdp_mul(f, g, u, O, K):
"""Multiply distributed polynomials in K[X]. """
if sdp_term_p(f):
if not f:
return f
else:
return sdp_mul_term(g, f[0], u, O, K)

if sdp_term_p(g):
if not g:
return g
else:
return sdp_mul_term(f, g[0], u, O, K)

h = {}

for fm, fc in f:
for gm, gc in g:
monom = monomial_mul(fm, gm)
coeff = fc * gc

if monom in h:
coeff += h[monom]

if not coeff:
del h[monom]
continue

h[monom] = coeff

return sdp_from_dict(h, O)

[docs]def sdp_sqr(f, u, O, K):
"""Square a distributed polynomial in K[X]. """
h = {}

for fm, fc in f:
for Fm, Fc in f:
monom = monomial_mul(fm, Fm)
coeff = fc * Fc

if monom in h:
coeff += h[monom]

if not coeff:
del h[monom]
continue

h[monom] = coeff

return sdp_from_dict(h, O)

[docs]def sdp_pow(f, n, u, O, K):
"""Raise f to the n-th power in K[X]. """
if not n:
return sdp_one(u, K)
if n < 0:
raise ValueError("can't raise a polynomial to negative power")
if n == 1 or not f or sdp_one_p(f, u, K):
return f

g = sdp_one(u, K)

while True:
n, m = n // 2, n

if m & 1:
g = sdp_mul(g, f, u, O, K)

if not n:
break

f = sdp_sqr(f, u, O, K)

return g

[docs]def sdp_monic(f, K):
"""Divides all coefficients by LC(f) in K[X]. """
if not f:
return f

lc_f = sdp_LC(f, K)

if K.is_one(lc_f):
return f
else:
return [ (m, K.quo(c, lc_f)) for m, c in f ]

[docs]def sdp_content(f, K):
"""Returns GCD of coefficients in K[X]. """
from sympy.polys.domains import QQ

if K.has_Field:
return K.one
else:
cont = K.zero

if K == QQ:
for c in f:
cont = K.gcd(cont, c)
else:
for c in f:
cont = K.gcd(cont, c)

if K.is_one(cont) and i:
break

return cont

[docs]def sdp_primitive(f, K):
"""Returns content and a primitive polynomial in K[X]. """
if K.has_Field:
return K.one, f
else:
cont = sdp_content(f, K)

if K.is_one(cont):
return cont, f
else:
return cont, [ (m, K.quo(c, cont)) for m, c in f ]

def _term_rr_div(a, b, K):
"""Division of two terms in over a ring. """
a_lm, a_lc = a
b_lm, b_lc = b

monom = monomial_div(a_lm, b_lm)

if not (monom is None or a_lc % b_lc):
return monom, K.quo(a_lc, b_lc)
else:
return None

def _term_ff_div(a, b, K):
"""Division of two terms in over a field. """
a_lm, a_lc = a
b_lm, b_lc = b

monom = monomial_div(a_lm, b_lm)

if monom is not None:
return monom, K.quo(a_lc, b_lc)
else:
return None

[docs]def sdp_div(f, G, u, O, K):
"""
Generalized polynomial division with remainder in K[X].

Given polynomial f and a set of polynomials g = (g_1, ..., g_n)
compute a set of quotients q = (q_1, ..., q_n) and remainder r
such that f = q_1*f_1 + ... + q_n*f_n + r, where r = 0 or r
is a completely reduced polynomial with respect to g.

References
==========

1. [Cox97]_
2. [Ajwa95]_

"""
Q, r = [ [] for _ in range(len(G)) ], []

if K.has_Field:
term_div = _term_ff_div
else:
term_div = _term_rr_div

while f:
for i, g in enumerate(G):
tq = term_div(sdp_LT(f, u, K), sdp_LT(g, u, K), K)

if tq is not None:
Q[i] = sdp_add_term(Q[i], tq, u, O, K)
f = sdp_sub(f, sdp_mul_term(g, tq, u, O, K), u, O, K)

break
else:
r = sdp_add_term(r, sdp_LT(f, u, K), u, O, K)
f = sdp_del_LT(f)

return Q, r

[docs]def sdp_rem(f, G, u, O, K):
"""Returns polynomial remainder in K[X]. """
r = {}

if K.has_Field:
term_div = _term_ff_div
else:
term_div = _term_rr_div

ltf = sdp_LT(f, u, K)
f = dict(f)
get = f.get
while f:
for g in G:
tq = term_div(ltf, sdp_LT(g, u, K), K)

if tq is not None:
m, c = tq
for mg, cg in g:
m1 = monomial_mul(mg, m)
c1 = get(m1, 0) - c*cg
if not c1:
del f[m1]
else:
f[m1] = c1
if f:
if O == lex:
ltm = max(f)
else:
ltm = max(f, key=lambda mx: O(mx))
ltf = ltm, f[ltm]

break
else:
ltm, ltc = ltf
if ltm in r:
r[ltm] += ltc
else:
r[ltm] = ltc
del f[ltm]
if f:
if O == lex:
ltm = max(f)
else:
ltm = max(f, key=lambda mx: O(mx))
ltf = ltm, f[ltm]
return sdp_from_dict(r, O)

[docs]def sdp_quo(f, g, u, O, K):
"""Returns polynomial quotient in K[x]. """
return sdp_div(f, g, u, O, K)[0]

[docs]def sdp_exquo(f, g, u, O, K):
"""Returns exact polynomial quotient in K[X]. """
q, r = sdp_div(f, g, u, O, K)

if not r:
return q
else:
raise ExactQuotientFailed(f, g)

[docs]def sdp_lcm(f, g, u, O, K):
"""
Computes LCM of two polynomials in K[X].

The LCM is computed as the unique generater of the intersection
of the two ideals generated by f and g. The approach is to
compute a Groebner basis with respect to lexicographic ordering
of t*f and (1 - t)*g, where t is an unrelated variable and
then filtering out the solution that doesn't contain t.

References
==========

1. [Cox97]_

"""
from sympy.polys.groebnertools import sdp_groebner

if not f or not g:
return []

if sdp_term_p(f) and sdp_term_p(g):
monom = monomial_lcm(sdp_LM(f, u), sdp_LM(g, u))

fc, gc = sdp_LC(f, K), sdp_LC(g, K)

if K.has_Field:
coeff = K.one
else:
coeff = K.lcm(fc, gc)

return [(monom, coeff)]

if not K.has_Field:
lcm = K.one
else:
fc, f = sdp_primitive(f, K)
gc, g = sdp_primitive(g, K)

lcm = K.lcm(fc, gc)

f_terms = tuple( ((1,) + m,  c) for m, c in f )
g_terms = tuple( ((0,) + m,  c) for m, c in g ) \
+ tuple( ((1,) + m, -c) for m, c in g )

F = sdp_sort(f_terms, lex)
G = sdp_sort(g_terms, lex)

basis = sdp_groebner([F, G], u, lex, K)

H = [ h for h in basis if sdp_indep_p(h, 0, u) ]

if K.is_one(lcm):
h = [ (m[1:], c)     for m, c in H[0] ]
else:
h = [ (m[1:], c * lcm) for m, c in H[0] ]

return sdp_sort(h, O)

[docs]def sdp_gcd(f, g, u, O, K):
"""Compute GCD of two polynomials in K[X] via LCM. """
if not K.has_Field:
fc, f = sdp_primitive(f, K)
gc, g = sdp_primitive(g, K)

gcd = K.gcd(fc, gc)

h = sdp_quo(sdp_mul(f, g, u, O, K),
sdp_lcm(f, g, u, O, K), u, O, K)

if not K.has_Field:
if K.is_one(gcd):
return h
else:
return [ (m, c * gcd) for m, c in h ]
else:
return sdp_monic(h, K)