Ticket #4539: extplural-more2.patch

sage/rings/polynomial/plural.pyx

@@ -20 +20 @@
 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 +111 @@
             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.<x,y,z> = 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.<x,y,z> = 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.<x,y,z> = 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  <Parent>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 +254 @@
            - 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.<x,y,z> = 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.<x,y,z> = 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.<x,y,z> = 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 (<Parent>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 (<ParentWithGens>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.<x,y,z> = 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.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex')
+##            sage: R == P
+##            True
+           
+##            sage: c[0,1] = -2
+##            sage: R.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex')
+##            sage: P == R
+##            False
+##         """
+###         return (<Parent>left)._richcmp(right, op)
+          ##return (<ParentWithGens>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.<x,y,z> = 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.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex')
+##            sage: R == P
+##            True
+           
+##            sage: c[0,1] = -2
+##            sage: R.<x,y,z> = 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 +582 @@
         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.<x,y,z> = 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.<x,y,z> = 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.<x,y,z> = 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.<x,y,z> = 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 <ParentWithBase>self is f._parent:
+            f = self._coerce_c(f)
+        if not <ParentWithBase>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.<x,y,z> = 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.<x,y,z> = 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.<x,y,z> = 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 = (<NCPolynomialRing_plural>a._parent)._coerce_c(b)
+
+        _a = a._poly
+        _b = b._poly
+        _r = (<NCPolynomialRing_plural>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.<x,y,z> = 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.<x,y,z> = 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.<x,y,z> = 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 <ParentWithBase>self is f._parent:
+            f = self._coerce_c(f)
+        if not <ParentWithBase>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.<x,y,z> = 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.<x,y,z> = 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.<x,y,z> = 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 (<NCPolynomial_plural>g) \
+                   and p_LmDivisibleBy((<NCPolynomial_plural>g)._poly, m, r):
+                flt = pDivide(f._poly, (<NCPolynomial_plural>g)._poly)
+                #p_SetCoeff(flt, n_Div( p_GetCoeff(f._poly, r) , p_GetCoeff((<NCPolynomial_plural>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.<x,y,z> = 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.<x,y,z> = QQ[]
+            sage: P.<x,y,z> = 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 = (<MPolynomialRing_libsingular>h._parent)._coerce_c(g)
+
+        r = (<MPolynomialRing_libsingular>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.<x,y,z> = 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 +1622 @@
     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