diff -r d806d30df184 sage/rings/polynomial/plural.pyx --- a/sage/rings/polynomial/plural.pyx Tue Jul 20 16:22:29 2010 +0200 +++ b/sage/rings/polynomial/plural.pyx Wed Jul 21 12:21:55 2010 +0200 @@ -20,10 +20,16 @@ from term_order import TermOrder -from sage.rings.polynomial.multi_polynomial_libsingular cimport MPolynomialRing_libsingular +from sage.rings.polynomial.multi_polynomial_libsingular cimport MPolynomialRing_libsingular +#from sage.rings.polynomial.multi_polynomial_libsingular cimport addwithcarry from sage.rings.polynomial.multi_polynomial_ring_generic import MPolynomialRing_generic +from sage.structure.parent cimport Parent +from sage.structure.element cimport CommutativeRingElement +from sage.rings.finite_rings.finite_field_prime_modn import FiniteField_prime_modn +from sage.rings.integer_ring import is_IntegerRing, ZZ + cdef class NCPolynomialRing_plural(Ring): def __init__(self, base_ring, n, names, c, d, order='degrevlex', check = True): order = TermOrder(order,n) @@ -105,18 +111,142 @@ sage: _ = gc.collect() """ singular_ring_delete(self._ring) + + def _element_constructor_(self, element): + """ + Make sure element is a valid member of self, and return the constructed element. + + EXAMPLES:: + sage: A. = FreeAlgebra(QQ, 3) + + sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex') + + We can construct elements from the base ring:: + + sage: P(1/2) + 1/2 + + + and all kinds of integers:: + + sage: P(17) + 17 + + sage: P(int(19)) + 19 + + sage: P(long(19)) + 19 + + TESTS:: + + Check conversion from self:: + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex') + sage: P.inject_variables() + Defining x, y, z + + sage: P(1/2) + 1/2 + + sage: P(x*y) + x*y + + sage: P(y*x) + -x*y + + Raw use of this class:: + sage: from sage.matrix.constructor import Matrix + sage: c = Matrix(3) + sage: c[0,1] = -2 + sage: c[0,2] = 1 + sage: c[1,2] = 1 + + sage: d = Matrix(3) + sage: d[0, 1] = 17 + + sage: from sage.rings.polynomial.plural import NCPolynomialRing_plural + sage: R. = NCPolynomialRing_plural(QQ, 3, c = c, d = d, order='lex') + sage: R(x*y) + x*y + + sage: P(17) + 17 + + sage: P(int(19)) + 19 + + Testing special cases:: + sage: P(1) + 1 + + sage: P(0) + 0 + """ - - - def _element_constructor_(self, x): - """ - Make sure x is a valid member of self, and return the constructed element. - """ - if x==0: - return self._zero_element - if x==1: - return self._one_element - raise NotImplementedError("not able to constructor "+repr(x)) #### ?????? + if element == 0: + return self._zero_element + if element == 1: + return self._one_element + + cdef poly *_p + cdef ring *_ring + cdef number *_n + + _ring = self._ring + + base_ring = self.base_ring() + + if(_ring != currRing): rChangeCurrRing(_ring) + + + if PY_TYPE_CHECK(element, CommutativeRingElement): + # base ring elements + if element.parent() is base_ring: + # shortcut for GF(p) + if isinstance(base_ring, FiniteField_prime_modn): + _p = p_ISet(int(element) % _ring.ch, _ring) + else: + _n = sa2si(element,_ring) + _p = p_NSet(_n, _ring) + + # also accepting ZZ + elif is_IntegerRing(element.parent()): + if isinstance(base_ring, FiniteField_prime_modn): + _p = p_ISet(int(element),_ring) + else: + _n = sa2si(base_ring(element),_ring) + _p = p_NSet(_n, _ring) + else: + # fall back to base ring + element = base_ring._coerce_c(element) + _n = sa2si(element,_ring) + _p = p_NSet(_n, _ring) + + # Accepting int + elif PY_TYPE_CHECK(element, int): + if isinstance(base_ring, FiniteField_prime_modn): + _p = p_ISet(int(element) % _ring.ch,_ring) + else: + _n = sa2si(base_ring(element),_ring) + _p = p_NSet(_n, _ring) + + # and longs + elif PY_TYPE_CHECK(element, long): + if isinstance(base_ring, FiniteField_prime_modn): + element = element % self.base_ring().characteristic() + _p = p_ISet(int(element),_ring) + else: + _n = sa2si(base_ring(element),_ring) + _p = p_NSet(_n, _ring) + + else: + raise NotImplementedError("not able to constructor "+repr(element) + + " of type "+ repr(type(element))) #### ?????? + + return new_NCP(self,_p) + + cpdef _coerce_map_from_(self, S): """ @@ -124,9 +254,182 @@ - the integer ring - other localizations away from fewer primes """ - if S is ZZ: + + if self.base_ring().has_coerce_map_from(S): return True + + + def __hash__(self): + """ + Return a hash for this noncommutative ring, that is, a hash of the string + representation of this polynomial ring. + + EXAMPLES:: + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex') + sage: hash(P) # somewhat random output + ... + + TESTS:: + + Check conversion from self:: + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex') + sage: from sage.matrix.constructor import Matrix + sage: c = Matrix(3) + sage: c[0,1] = -1 + sage: c[0,2] = 1 + sage: c[1,2] = 1 + + sage: from sage.rings.polynomial.plural import NCPolynomialRing_plural + sage: R. = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex') + sage: hash(R) == hash(P) + True + """ + return hash(self.__repr__()) + +## def _richcmp_(self, right, int op): +## return (self)._richcmp_helper(right, op) + + +## ## cdef int _cmp_c_impl(left, Parent right) except -2: +## ## return 0#cmp(type(left),type(right)) + +## def __cmp__(self, x): +## """ +## EXAMPLES:: + +## """ +## if PY_TYPE_CHECK(x, NCPolynomialRing_plural): +## print "huhu" +## return 0 +## print "hihi" +## return 0# return cmp(type(self), type(x)) + +## def __richcmp__(self, right, int op): +## print "CMP" +## print "base", self.base_ring() +## print "gens" +## for elt in self.gens(): +## print elt +## print self.term_order() +## for elt in self.relations(): +## print elt +## print PY_TYPE_CHECK(right, NCPolynomialRing_plural), "huhu" + +## print "base", right.base_ring() +## print "gens" +## for elt in right.gens(): +## print elt +## print right.term_order() +## for elt in right.relations(): +## print elt +## return False + +## print cmp(self.base_ring(), right.base_ring()) +## print cmp( map(str, self.gens()), map(str, right.gens())) +## print cmp( map(str, self.relations()), map(str, right.relations())) + +## return False +## if PY_TYPE_CHECK(right, NCPolynomialRing_plural): +## return cmp( (self.base_ring(), map(str, self.gens()), +## self.term_order(), map(str, self.relations())), +## (right.base_ring(), map(str, right.gens()), +## right.term_order(), map(str, right.relations())) +## ) +## else: +## return cmp(type(self),type(right)) + +## # return False + + # def __richcmp__(left, right, int op): + # return (left)._richcmp(right, op) + # return False +## r""" +## Multivariate polynomial rings are said to be equal if: + +## - their base rings match, +## - their generator names match, +## - their term orderings match, and +## - their relations match. + + +## EXAMPLES:: +## sage: A. = FreeAlgebra(QQ, 3) +## sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex') + +## sage: P == P +## True +## sage: Q = copy(P) +## sage: Q == P +## True + +## sage: from sage.matrix.constructor import Matrix +## sage: c = Matrix(3) +## sage: c[0,1] = -1 +## sage: c[0,2] = 1 +## sage: c[1,2] = 1 +## sage: from sage.rings.polynomial.plural import NCPolynomialRing_plural +## sage: R. = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex') +## sage: R == P +## True + +## sage: c[0,1] = -2 +## sage: R. = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex') +## sage: P == R +## False +## """ +### return (left)._richcmp(right, op) + ##return (left)._richcmp(right, op) +## cdef int _cmp_c_impl(left, Parent right) except -2: +## r""" +## Multivariate polynomial rings are said to be equal if: + +## - their base rings match, +## - their generator names match, +## - their term orderings match, and +## - their relations match. + + +## EXAMPLES:: +## sage: A. = FreeAlgebra(QQ, 3) +## sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex') + +## sage: P == P +## True +## sage: Q = copy(P) +## sage: Q == P +## True + +## sage: from sage.matrix.constructor import Matrix +## sage: c = Matrix(3) +## sage: c[0,1] = -1 +## sage: c[0,2] = 1 +## sage: c[1,2] = 1 +## sage: from sage.rings.polynomial.plural import NCPolynomialRing_plural +## sage: R. = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex') +## sage: R == P +## True + +## sage: c[0,1] = -2 +## sage: R. = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex') +## sage: P == R +## False +## """ +## print "huhu", PY_TYPE_CHECK(right, NCPolynomialRing_plural) +## if PY_TYPE_CHECK(right, NCPolynomialRing_plural): +## return True +## return cmp( (left.base_ring(), map(str, left.gens()), +## left.term_order(), map(str, left.relations())), +## (right.base_ring(), map(str, right.gens()), +## right.term_order(), map(str, right.relations())) +## ) +## else: +## return cmp(type(left),type(right)) + + + def __pow__(self, n, _): """ Return the free module of rank `n` over this ring. @@ -279,6 +582,405 @@ W = self._list_to_ring(L) return new_NRing(W, self.base_ring()) + + ### The following methods are handy for implementing Groebner + ### basis algorithms. They do only superficial type/sanity checks + ### and should be called carefully. + + def monomial_quotient(self, NCPolynomial_plural f, NCPolynomial_plural g, coeff=False): + r""" + Return ``f/g``, where both ``f`` and`` ``g`` are treated as + monomials. + + Coefficients are ignored by default. + + INPUT: + + - ``f`` - monomial + - ``g`` - monomial + - ``coeff`` - divide coefficients as well (default: ``False``) + + EXAMPLES:: + + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + + sage: P.monomial_quotient(3/2*x*y,x,coeff=True) + 3/2*y + + Note, that `\ZZ` behaves different if ``coeff=True``:: + + sage: P.monomial_quotient(2*x,3*x) + 1 + + sage: A. = FreeAlgebra(ZZ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + + sage: P.monomial_quotient(2*x,3*x,coeff=True) + Traceback (most recent call last): + ... + ArithmeticError: Cannot divide these coefficients. + + TESTS:: + sage: A. = FreeAlgebra(QQ, 3) + sage: R = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: R.inject_variables() + Defining x, y, z + + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + + sage: P.monomial_quotient(x*y,x) + y + +## sage: P.monomial_quotient(x*y,R.gen()) +## y + + sage: P.monomial_quotient(P(0),P(1)) + 0 + + sage: P.monomial_quotient(P(1),P(0)) + Traceback (most recent call last): + ... + ZeroDivisionError + + sage: P.monomial_quotient(P(3/2),P(2/3), coeff=True) + 9/4 + + sage: P.monomial_quotient(x,y) # Note the wrong result + x*y^1048575*z^1048575 # 64-bit + x*y^65535*z^65535 # 32-bit + + sage: P.monomial_quotient(x,P(1)) + x + + .. warning:: + + Assumes that the head term of f is a multiple of the head + term of g and return the multiplicant m. If this rule is + violated, funny things may happen. + """ + cdef poly *res + cdef ring *r = self._ring + cdef number *n, *denom + + if not self is f._parent: + f = self._coerce_c(f) + if not self is g._parent: + g = self._coerce_c(g) + + if(r != currRing): rChangeCurrRing(r) + + if not f._poly: + return self._zero_element + if not g._poly: + raise ZeroDivisionError + + res = pDivide(f._poly,g._poly) + if coeff: + if r.ringtype == 0 or r.cf.nDivBy(p_GetCoeff(f._poly, r), p_GetCoeff(g._poly, r)): + n = r.cf.nDiv( p_GetCoeff(f._poly, r) , p_GetCoeff(g._poly, r)) + p_SetCoeff0(res, n, r) + else: + raise ArithmeticError("Cannot divide these coefficients.") + else: + p_SetCoeff0(res, n_Init(1, r), r) + return new_NCP(self, res) + + def monomial_divides(self, NCPolynomial_plural a, NCPolynomial_plural b): + """ + Return ``False`` if a does not divide b and ``True`` + otherwise. + + Coefficients are ignored. + + INPUT: + + - ``a`` -- monomial + + - ``b`` -- monomial + + EXAMPLES:: + + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + + sage: P.monomial_divides(x*y*z, x^3*y^2*z^4) + True + sage: P.monomial_divides(x^3*y^2*z^4, x*y*z) + False + + TESTS:: + sage: A. = FreeAlgebra(QQ, 3) + sage: Q = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: Q.inject_variables() + Defining x, y, z + + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + + sage: P.monomial_divides(P(1), P(0)) + True + sage: P.monomial_divides(P(1), x) + True + """ + cdef poly *_a + cdef poly *_b + cdef ring *_r + if a._parent is not b._parent: + b = (a._parent)._coerce_c(b) + + _a = a._poly + _b = b._poly + _r = (a._parent)._ring + + if _a == NULL: + raise ZeroDivisionError + if _b == NULL: + return True + + if not p_DivisibleBy(_a, _b, _r): + return False + else: + return True + + + def monomial_lcm(self, NCPolynomial_plural f, NCPolynomial_plural g): + """ + LCM for monomials. Coefficients are ignored. + + INPUT: + + - ``f`` - monomial + + - ``g`` - monomial + + EXAMPLES:: + + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + + sage: P.monomial_lcm(3/2*x*y,x) + x*y + + TESTS:: + + sage: A. = FreeAlgebra(QQ, 3) + sage: R = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: R.inject_variables() + Defining x, y, z + + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + +## sage: P.monomial_lcm(x*y,R.gen()) +## x*y + + sage: P.monomial_lcm(P(3/2),P(2/3)) + 1 + + sage: P.monomial_lcm(x,P(1)) + x + """ + cdef poly *m = p_ISet(1,self._ring) + + if not self is f._parent: + f = self._coerce_c(f) + if not self is g._parent: + g = self._coerce_c(g) + + if f._poly == NULL: + if g._poly == NULL: + return self._zero_element + else: + raise ArithmeticError, "Cannot compute LCM of zero and nonzero element." + if g._poly == NULL: + raise ArithmeticError, "Cannot compute LCM of zero and nonzero element." + + if(self._ring != currRing): rChangeCurrRing(self._ring) + + pLcm(f._poly, g._poly, m) + p_Setm(m, self._ring) + return new_NCP(self,m) + + def monomial_reduce(self, NCPolynomial_plural f, G): + """ + Try to find a ``g`` in ``G`` where ``g.lm()`` divides + ``f``. If found ``(flt,g)`` is returned, ``(0,0)`` otherwise, + where ``flt`` is ``f/g.lm()``. + + It is assumed that ``G`` is iterable and contains *only* + elements in this polynomial ring. + + Coefficients are ignored. + + INPUT: + + - ``f`` - monomial + - ``G`` - list/set of mpolynomials + + EXAMPLES:: + + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + + sage: f = x*y^2 + sage: G = [ 3/2*x^3 + y^2 + 1/2, 1/4*x*y + 2/7, 1/2 ] + sage: P.monomial_reduce(f,G) + (y, 1/4*x*y + 2/7) + + TESTS:: + sage: A. = FreeAlgebra(QQ, 3) + sage: Q = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: Q.inject_variables() + Defining x, y, z + + sage: A. = FreeAlgebra(QQ, 3) + sage: P = A.g_algebra(relations={y*x:-x*y}, order='lex') + sage: P.inject_variables() + Defining x, y, z + sage: f = x*y^2 + sage: G = [ 3/2*x^3 + y^2 + 1/2, 1/4*x*y + 2/7, 1/2 ] + + sage: P.monomial_reduce(P(0),G) + (0, 0) + + sage: P.monomial_reduce(f,[P(0)]) + (0, 0) + """ + cdef poly *m = f._poly + cdef ring *r = self._ring + cdef poly *flt + + if not m: + return f,f + + for g in G: + if PY_TYPE_CHECK(g, NCPolynomial_plural) \ + and (g) \ + and p_LmDivisibleBy((g)._poly, m, r): + flt = pDivide(f._poly, (g)._poly) + #p_SetCoeff(flt, n_Div( p_GetCoeff(f._poly, r) , p_GetCoeff((g)._poly, r), r), r) + p_SetCoeff(flt, n_Init(1, r), r) + return new_NCP(self,flt), g + return self._zero_element,self._zero_element + + def monomial_pairwise_prime(self, NCPolynomial_plural g, NCPolynomial_plural h): + """ + Return ``True`` if ``h`` and ``g`` are pairwise prime. Both + are treated as monomials. + + Coefficients are ignored. + + INPUT: + + - ``h`` - monomial + - ``g`` - monomial + + EXAMPLES:: + + sage: P. = QQ[] + sage: P.monomial_pairwise_prime(x^2*z^3, y^4) + True + + sage: P.monomial_pairwise_prime(1/2*x^3*y^2, 3/4*y^3) + False + + TESTS:: + + sage: Q. = QQ[] + sage: P. = QQ[] + sage: P.monomial_pairwise_prime(x^2*z^3, Q('y^4')) + True + + sage: P.monomial_pairwise_prime(1/2*x^3*y^2, Q(0)) + True + + sage: P.monomial_pairwise_prime(P(1/2),x) + False + """ + cdef int i + cdef ring *r + cdef poly *p, *q + + if h._parent is not g._parent: + g = (h._parent)._coerce_c(g) + + r = (h._parent)._ring + p = g._poly + q = h._poly + + if p == NULL: + if q == NULL: + return False #GCD(0,0) = 0 + else: + return True #GCD(x,0) = 1 + + elif q == NULL: + return True # GCD(0,x) = 1 + + elif p_IsConstant(p,r) or p_IsConstant(q,r): # assuming a base field + return False + + for i from 1 <= i <= r.N: + if p_GetExp(p,i,r) and p_GetExp(q,i,r): + return False + return True + + def monomial_all_divisors(self, NCPolynomial_plural t): + """ + Return a list of all monomials that divide ``t``. + + Coefficients are ignored. + + INPUT: + + - ``t`` - a monomial + + OUTPUT: + a list of monomials + + + EXAMPLES:: + + sage: P. = QQ[] + sage: P.monomial_all_divisors(x^2*z^3) + [x, x^2, z, x*z, x^2*z, z^2, x*z^2, x^2*z^2, z^3, x*z^3, x^2*z^3] + + ALGORITHM: addwithcarry idea by Toon Segers + """ + + M = list() + + cdef ring *_ring = self._ring + cdef poly *maxvector = t._poly + cdef poly *tempvector = p_ISet(1, _ring) + + pos = 1 + + while not p_ExpVectorEqual(tempvector, maxvector, _ring): + tempvector = addwithcarry(tempvector, maxvector, pos, _ring) + M.append(new_NCP(self, p_Copy(tempvector,_ring))) + return M + + + cdef class NCPolynomial_plural(RingElement): def __init__(self, NCPolynomialRing_plural parent): self._poly = NULL @@ -920,3 +1622,11 @@ I = H.ideal([H.gen(i) *H.gen(i) for i in alt_vars]).twostd() return H.quotient(I) +cdef poly *addwithcarry(poly *tempvector, poly *maxvector, int pos, ring *_ring): + if p_GetExp(tempvector, pos, _ring) < p_GetExp(maxvector, pos, _ring): + p_SetExp(tempvector, pos, p_GetExp(tempvector, pos, _ring)+1, _ring) + else: + p_SetExp(tempvector, pos, 0, _ring) + tempvector = addwithcarry(tempvector, maxvector, pos + 1, _ring) + p_Setm(tempvector, _ring) + return tempvector