Add FGLM #538

Merged
merged 3 commits into from Aug 12, 2011
View
2 sympy/polys/__init__.py
@@ -25,7 +25,7 @@
real_roots, nroots, ground_roots,
nth_power_roots_poly,
cancel,
- reduced, groebner,
+ reduced, groebner, fglm
)
from polyfuncs import (
View
196 sympy/polys/groebnertools.py
@@ -12,17 +12,10 @@
)
from sympy.polys.distributedpolys import (
- sdp_LC,
- sdp_LM,
- sdp_LT,
- sdp_mul_term,
- sdp_sub,
- sdp_mul_term,
- sdp_monic,
- sdp_rem,
- sdp_strip,
- _term_ff_div,
- _term_rr_div,
+ sdp_LC, sdp_LM, sdp_LT, sdp_mul_term,
+ sdp_sub, sdp_mul_term, sdp_monic,
+ sdp_rem, sdp_strip, sdp_sort,
+ _term_ff_div, _term_rr_div,
)
from sympy.polys.polyerrors import (
@@ -827,3 +820,184 @@ def monomial_divides(m1, m2):
return False
return True
+
+# FGLM
+
+
+def is_zero_dimensional(F, u, O, K):
+ """
+ Checks if the ideal generated by ``F`` is zero-dimensional.
+
+ The algorithm checks if the set of monomials not divisible by a
+ leading monomial of any element of ``F`` is bounded. In general
+ ``F`` has to be a Groebner basis w.r.t. ``O`` but if ``True``
+ is returned, then the ideal is zero-dimensional.
+
+ **References**
+ David A. Cox, John B. Little, Donal O'Shea. Ideals, Varieties
+ and Algorithms, 3rd edition, p. 234
+ """
+ def single_var(m):
+ n = 0
+ for e in m:
+ if e != 0:
+ n += 1
+ return n == 1
+
+ exponents = (0,) * (u + 1)
+
+ for f in F:
+ if single_var(sdp_LM(f, u)):
+ exponents = monomial_mul(exponents, sdp_LM(f, O)) # == sum of exponent vectors
+
+ product = 1
+ for e in exponents:
+ product *= e
+
+ # If product == 0, then there's a variable for which there's
+ # no degree bound.
+ return product != 0
+
+
+def matrix_fglm(F, u, O_from, O_to, K):
+ """
+ Converts the reduced Groebner basis ``F`` of a zero-dimensional
+ ideal w.r.t. ``O_from`` to a reduced Groebner basis
+ w.r.t. ``O_to``.
+
+ **References**
+ J.C. Faugere, P. Gianni, D. Lazard, T. Mora (1994). Efficient
+ Computation of Zero-dimensional Groebner Bases by Change of
+ Ordering
+
+ J.C. Faugere's lecture notes:
+ http://www-salsa.lip6.fr/~jcf/Papers/2010_MPRI5e.pdf
+ """
+ old_basis = _basis(F, u, O_from, K)
+ M = _representing_matrices(old_basis, F, u, O_from, K)
+
+ # V contains the normalforms (wrt O_from) of S
+ S = [(0,) * (u + 1)]
+ V = [[K.one] + [K.zero] * (len(old_basis) - 1)]
+ G = []
+
+ L = [(i, 0) for i in xrange(u + 1)] # (i, j) corresponds to x_i * S[j]
+ L.sort(key=lambda (k, l): O_to(_incr_k(S[l], k)), reverse=True)
+ t = L.pop()
+
+ P = _identity_matrix(len(old_basis), K)
+
+ while True:
+ s = len(S)
+ v = _matrix_mul(M[t[0]], V[t[1]], K)
+ _lambda = _matrix_mul(P, v, K)
+
+ if all([_lambda[i] == K.zero for i in xrange(s, len(old_basis))]):
+ # there is a linear combination of v by V
+
+ lt = [(_incr_k(S[t[1]], t[0]), K.one)]
+ rest = sdp_strip(sdp_sort([(S[i], _lambda[i]) for i in xrange(s)], O_to))
+ g = sdp_sub(lt, rest, u, O_to, K)
+
+ if g != []:
+ G.append(g)
+
+ else:
+ # v is linearly independant from V
+ P = _update(s, _lambda, P, K)
+ S.append(_incr_k(S[t[1]], t[0]))
+ V.append(v)
+
+ L.extend([(i, s) for i in xrange(u + 1)])
+ L = list(set(L))
+ L.sort(key=lambda (k, l): O_to(_incr_k(S[l], k)), reverse=True)
+
+ L = [(k, l) for (k, l) in L if \
+ all([monomial_div(_incr_k(S[l], k), sdp_LM(g, u)) is None for g in G])]
+
+ if L == []:
+ return sorted(G, key=lambda g: O_to(sdp_LM(g, u)), reverse=True)
+
+ t = L.pop()
+
+
+def _incr_k(m, k):
+ return tuple(list(m[:k]) + [m[k] + 1] + list(m[k + 1:]))
+
+
+def _identity_matrix(n, K):
+ M = [[K.zero] * n for _ in xrange(n)]
+
+ for i in xrange(n):
+ M[i][i] = K.one
+
+ return M
+
+
+def _matrix_mul(M, v, K):
+ return [sum([row[i] * v[i] for i in xrange(len(v))]) for row in M]
+
+
+def _update(s, _lambda, P, K):
+ """
+ Update ``P`` such that for the updated `P'` `P' v = e_{s}`.
+ """
+ k = min([j for j in xrange(s, len(_lambda)) if _lambda[j] != 0])
+
+ for r in xrange(len(_lambda)):
+ if r != k:
+ P[r] = [P[r][j] - (P[k][j] * _lambda[r]) / _lambda[k] for j in xrange(len(P[r]))]
+
+ P[k] = [P[k][j] / _lambda[k] for j in xrange(len(P[k]))]
+
+ P[k], P[s] = P[s], P[k]
+
+ return P
+
+
+def _representing_matrices(basis, G, u, O, K):
+ """
+ Compute the matrices corresponding to the linear maps `m \mapsto
+ x_i m` for all variables `x_i`.
+ """
+ def var(i):
+ return tuple([0] * i + [1] + [0] * (u - i))
+
+ def representing_matrix(m):
+ M = [[K.zero] * len(basis) for _ in xrange(len(basis))]
+
+ for i, v in enumerate(basis):
+ r = sdp_rem([(monomial_mul(m, v), K.one)], G, u, O, K)
+
+ for term in r:
+ j = basis.index(term[0])
+ M[j][i] = term[1]
+
+ return M
+
+ return [representing_matrix(var(i)) for i in xrange(u + 1)]
+
+
+def _basis(G, u, O, K):
+ """
+ Computes a list of monomials which are not divisible by the leading
+ monomials wrt to ``O`` of ``G``. These monomials are a basis of
+ `K[X_1, \ldots, X_n]/(G)`.
+ """
+ leading_monomials = [sdp_LM(g, u) for g in G]
+ candidates = [(0,) * (u + 1)]
+ basis = []
+
+ while candidates:
+ t = candidates.pop()
+ basis.append(t)
+
+ new_candidates = [_incr_k(t, k) for k in xrange(u + 1) \
+ if all([monomial_div(_incr_k(t, k), lmg) is None \
+ for lmg in leading_monomials])]
+ candidates.extend(new_candidates)
+ candidates.sort(key=lambda m: O(m), reverse=True)
+
+ basis = list(set(basis))
+
+ return sorted(basis, key=lambda m: O(m))
View
74 sympy/polys/polytools.py
@@ -38,7 +38,8 @@
)
from sympy.polys.groebnertools import (
- sdp_groebner,
+ sdp_groebner, matrix_fglm,
+ is_zero_dimensional,
)
from sympy.polys.monomialtools import (
@@ -5345,6 +5346,77 @@ def groebner(F, *gens, **args):
else:
return G
+def fglm(G, order, *gens, **args):
+ """
+ Convert reduced Groebner basis ``G`` of zero-dimensional ideal from
+ ``from_order`` to ``to_order``.
+ ``fglm`` does not check if ``G`` is a Groebner basis or if
+ it is reduced but does check if the ideal generated by ``G``
+ is zero-dimensional.
+
+ **Examples**
+
+ >>> from sympy.abc import x, y
+ >>> from sympy import groebner, fglm
+ >>> F = [x**2 - 3*y - x + 1, y**2 - 2*x + y - 1]
+ >>> G = groebner(F, x, y, order='grlex')
+ >>> fglm(G, ('grlex', 'lex'), x, y)
+ [2*x - y**2 - y + 1, y**4 + 2*y**3 - 3*y**2 - 16*y + 7]
+ >>> groebner(F, x, y, order='lex')
+ [2*x - y**2 - y + 1, y**4 + 2*y**3 - 3*y**2 - 16*y + 7]
+
+ **References**
+ J.C. Faugere, P. Gianni, D. Lazard, T. Mora (1994). Efficient
+ Computation of Zero-dimensional Groebner Bases by Change of
+ Ordering
+
+ J.C. Faugere's lecture notes:
+ http://www-salsa.lip6.fr/~jcf/Papers/2010_MPRI5e.pdf
+ """
+ options.allowed_flags(args, ['polys'])
+
+ if type(order) == str:
+ from_order, to_order = order, 'lex'
+ elif type(order) == tuple and len(order) == 2:
+ from_order, to_order = order
+ elif type(order) == tuple and len(order) == 1:
+ from_order, to_order = order[0], 'lex'
+ else:
+ raise TypeError("order has to be a string or a tuple of at most two strings")
+
+ try:
+ polys, opt = parallel_poly_from_expr(G, *gens, **args)
+ except PolificationFailed, exc:
+ raise ComputationFailed('fglm', len(G), exc)
+
+ domain = opt.domain
+
+ if domain.has_assoc_Field:
+ opt.domain = domain.get_field()
+ else:
+ raise DomainError("can't convert Groebner basis over %s" % domain)
+
+ for i, poly in enumerate(polys):
+ poly = poly.set_domain(opt.domain).rep.to_dict()
+ polys[i] = sdp_from_dict(poly, monomial_key(from_order))
+
+ level = len(flatten(gens)) - 1
+
+ if not is_zero_dimensional(polys, level, monomial_key(from_order), opt.domain):
+ raise NotImplementedError("Can't convert Groebner bases of ideals with positive dimension")
+
+ G = matrix_fglm(polys, level, monomial_key(from_order), monomial_key(to_order), opt.domain)
+ opt.order = monomial_key(to_order)
+ G = [Poly._from_dict(dict(g), opt).primitive()[1] for g in G]
+
+ if not domain.has_Field:
+ G = [g.clear_denoms(convert=True)[1] for g in G]
+
+ if not opt.polys:
+ return [g.as_expr() for g in G]
+ else:
+ return G
+
def poly(expr, *gens, **args):
"""
Efficiently transform an expression into a polynomial.
View
52 sympy/polys/tests/test_groebnertools.py
@@ -5,22 +5,12 @@
)
from sympy.polys.groebnertools import (
- sdp_groebner,
- sig,
- sig_key,
- sig_cmp,
- lbp,
- lbp_cmp,
- lbp_key,
- critical_pair,
- cp_cmp,
- cp_key,
- is_rewritable_or_comparable,
- Sign,
- Polyn,
- Num,
- s_poly,
- f5_reduce,
+ sdp_groebner, sig, sig_key, sig_cmp,
+ lbp, lbp_cmp, lbp_key, critical_pair,
+ cp_cmp, cp_key, is_rewritable_or_comparable,
+ Sign, Polyn, Num, s_poly, f5_reduce,
+ matrix_fglm, is_zero_dimensional,
+ _basis, _representing_matrices,
)
from sympy.polys.monomialtools import (
@@ -438,3 +428,33 @@ def test_f5_reduce():
s = lbp(sig(Sign(s)[0], 100), Polyn(s), Num(s))
assert f5_reduce(s, F, 2, O_lex, QQ) == s
+def test_matrix_fglm():
+ pass # see test_polytools.py
+
+def test_is_zero_dimensional():
+ F = [[((3, 0), QQ.one), ((0, 2), QQ.one)]]
+
+ assert is_zero_dimensional(F, 1, O_lex, QQ) == False
+
+ F = [[((1, 0), QQ.one)], [((0, 1), QQ.one)]]
+
+ assert is_zero_dimensional(F, 1, O_lex, QQ) == True
+
+ F = [[((1, 0, 0, 0), QQ.one)], [((0, 1, 0, 0), QQ.one)], [((0, 0, 0, 1), QQ.one)]]
+
+ assert is_zero_dimensional(F, 3, O_grevlex, QQ) == False
+
+def test_representing_matrices():
+ basis = [(0, 0), (0, 1), (1, 0), (1, 1)]
+ F = [[((2, 0), QQ(1,1)), ((1, 0), QQ(-1,1)), ((0, 1), QQ(-3,1)), ((0, 0), QQ(1,1))],
+ [((0, 2), QQ(1,1)), ((1, 0), QQ(-2,1)), ((0, 1), QQ(1,1)), ((0, 0), QQ(-1,1))]]
+
+ assert _representing_matrices(basis, F, 1, O_grlex, QQ) ==[ \
+ [[QQ(0,1), QQ(0,1), QQ(-1,1), QQ(3,1)],
+ [QQ(0,1), QQ(0,1), QQ(3,1), QQ(-4,1)],
+ [QQ(1,1), QQ(0,1), QQ(1,1), QQ(6,1)],
+ [QQ(0,1), QQ(1,1), QQ(0,1), QQ(1,1)]],
+ [[QQ(0,1), QQ(1,1), QQ(0,1), QQ(-2,1)],
+ [QQ(1,1), QQ(-1,1), QQ(0,1), QQ(6,1)],
+ [QQ(0,1), QQ(2,1), QQ(0,1), QQ(3,1)],
+ [QQ(0,1), QQ(0,1), QQ(1,1), QQ(-1,1)]]]
View
29 sympy/polys/tests/test_polytools.py
@@ -24,7 +24,7 @@
real_roots, nroots, ground_roots,
nth_power_roots_poly,
cancel,
- reduced, groebner)
+ reduced, groebner, fglm)
from sympy.polys.polyerrors import (
MultivariatePolynomialError,
@@ -2412,6 +2412,33 @@ def test_groebner():
raises(DomainError, "groebner([x**2 + 2.0*y], x, y)")
raises(ComputationFailed, "groebner([1])")
+def test_fglm():
+ a, b, c, d = symbols('a b c d')
+ F = [a+b+c+d, a*b + a*d + b*c + b*d, a*b*c + a*b*d + a*c*d + b*c*d, a*b*c*d - 1]
+ G = groebner(F, a, b, c, d, order='grlex')
+
+ assert fglm(G, ('grlex'), a, b, c, d) == \
+ [4*a + 3*d**9 - 4*d**5 - 3*d, 4*b + 4*c - 3*d**9 + 4*d**5 + 7*d,
+ 4*c**2 + 3*d**10 - 4*d**6 - 3*d**2, 4*c*d**4 + 4*c - d**9 + 4*d**5 + 5*d,
+ d**12 - d**8 - d**4 + 1]
+
+ x, t = symbols('x t')
+ F = [9*x**8 + 36*x**7 - 32*x**6 - 252*x**5 - 78*x**4 + 468*x**3 + 288*x**2 - 108*x + 9,
+ -72*t*x**7 - 252*t*x**6 + 192*t*x**5 + 1260*t*x**4 + 312*t*x**3 - 404*t*x**2 - 576*t*x + \
+ 108*t - 72*x**7 - 256*x**6 + 192*x**5 + 1280*x**4 + 312*x**3 - 576*x + 96]
+ G = groebner(F, t, x, order='grlex')
+
+ assert fglm(G, 'grlex', t, x) == [203577793572507451707*t + 627982239411707112*x**7 - \
+ 666924143779443762*x**6 - 10874593056632447619*x**5 + 5119998792707079562*x**4 + \
+ 72917161949456066376*x**3 + 20362663855832380362*x**2 - 142079311455258371571*x + 183756699868981873194,
+ 9*x**8 + 36*x**7 - 32*x**6 - 252*x**5 - 78*x**4 + 468*x**3 + 288*x**2 - 108*x + 9]
+
+ x, y = symbols('x y')
+ F = [x**2 - x - 3*y + 1, -2*x + y**2 + y - 1]
+ G = groebner(F, x, y, order='lex')
+
+ assert fglm(G, ('lex', 'grlex'), x, y) == [x**2 - x - 3*y + 1, -2*x + y**2 + y - 1]
+
def test_poly():
assert poly(x) == Poly(x, x)
assert poly(y) == Poly(y, y)