# Ticket #4539: extplural-more2.patch

File extplural-more2.patch, 24.5 KB (added by AlexanderDreyer, 3 years ago)

noncommunative ring functionality

• ## sage/rings/polynomial/plural.pyx

`diff -r d806d30df184 sage/rings/polynomial/plural.pyx`
 a 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) 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): """ - 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. 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 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