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)