Ticket #4539: extplural-more2.patch

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

noncommunative ring functionality

  • sage/rings/polynomial/plural.pyx

    diff -r d806d30df184 sage/rings/polynomial/plural.pyx
    a b  
    2020from term_order import TermOrder
    2121
    2222
    23 from sage.rings.polynomial.multi_polynomial_libsingular cimport MPolynomialRing_libsingular
     23from sage.rings.polynomial.multi_polynomial_libsingular cimport MPolynomialRing_libsingular
     24#from sage.rings.polynomial.multi_polynomial_libsingular cimport addwithcarry
    2425from sage.rings.polynomial.multi_polynomial_ring_generic import MPolynomialRing_generic
    2526
    2627
     28from sage.structure.parent cimport Parent
     29from sage.structure.element cimport CommutativeRingElement
     30from sage.rings.finite_rings.finite_field_prime_modn import FiniteField_prime_modn
     31from sage.rings.integer_ring import is_IntegerRing, ZZ
     32
    2733cdef class NCPolynomialRing_plural(Ring):
    2834    def __init__(self, base_ring, n, names, c, d, order='degrevlex', check = True):
    2935        order = TermOrder(order,n)
     
    105111            sage: _ = gc.collect()
    106112        """
    107113        singular_ring_delete(self._ring)
     114   
     115    def _element_constructor_(self, element):
     116        """
     117        Make sure element is a valid member of self, and return the constructed element.
     118
     119        EXAMPLES::
     120            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     121
     122            sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex')
     123
     124        We can construct elements from the base ring::
     125
     126            sage: P(1/2)
     127            1/2
     128           
     129
     130        and all kinds of integers::
     131
     132            sage: P(17)
     133            17
     134
     135            sage: P(int(19))
     136            19
     137
     138            sage: P(long(19))
     139            19
     140           
     141        TESTS::
     142
     143        Check conversion from self::
     144            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     145            sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex')
     146            sage: P.inject_variables()
     147            Defining x, y, z
     148
     149            sage: P(1/2)
     150            1/2
     151
     152            sage: P(x*y)
     153            x*y
     154
     155            sage: P(y*x)
     156            -x*y         
     157
     158        Raw use of this class::
     159            sage: from sage.matrix.constructor  import Matrix
     160            sage: c = Matrix(3)
     161            sage: c[0,1] = -2
     162            sage: c[0,2] = 1
     163            sage: c[1,2] = 1
     164
     165            sage: d = Matrix(3)
     166            sage: d[0, 1] = 17
     167
     168            sage: from sage.rings.polynomial.plural import NCPolynomialRing_plural
     169            sage: R.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = d, order='lex')
     170            sage: R(x*y)
     171            x*y
     172           
     173            sage: P(17)
     174            17
     175
     176            sage: P(int(19))
     177            19
     178
     179        Testing special cases::
     180            sage: P(1)
     181            1
     182
     183            sage: P(0)
     184            0
     185        """
    108186       
    109        
    110    
    111     def _element_constructor_(self, x):
    112        """
    113        Make sure x is a valid member of self, and return the constructed element.
    114        """
    115        if x==0:
    116            return self._zero_element
    117        if x==1:
    118            return self._one_element
    119        raise NotImplementedError("not able to constructor "+repr(x)) #### ??????
     187        if element == 0:
     188            return self._zero_element
     189        if element == 1:
     190            return self._one_element
     191
     192        cdef poly *_p
     193        cdef ring *_ring
     194        cdef number *_n
     195       
     196        _ring = self._ring
     197       
     198        base_ring = self.base_ring()
     199       
     200        if(_ring != currRing): rChangeCurrRing(_ring)
     201
     202
     203        if PY_TYPE_CHECK(element, CommutativeRingElement):
     204            # base ring elements
     205            if  <Parent>element.parent() is base_ring:
     206                # shortcut for GF(p)
     207                if isinstance(base_ring, FiniteField_prime_modn):
     208                    _p = p_ISet(int(element) % _ring.ch, _ring)
     209                else:
     210                    _n = sa2si(element,_ring)
     211                    _p = p_NSet(_n, _ring)
     212                   
     213            # also accepting ZZ
     214            elif is_IntegerRing(element.parent()):
     215                if isinstance(base_ring, FiniteField_prime_modn):
     216                    _p = p_ISet(int(element),_ring)
     217                else:
     218                    _n = sa2si(base_ring(element),_ring)
     219                    _p = p_NSet(_n, _ring)
     220            else:
     221                # fall back to base ring
     222                element = base_ring._coerce_c(element)
     223                _n = sa2si(element,_ring)
     224                _p = p_NSet(_n, _ring)
     225
     226        # Accepting int
     227        elif PY_TYPE_CHECK(element, int):
     228            if isinstance(base_ring, FiniteField_prime_modn):
     229                _p = p_ISet(int(element) % _ring.ch,_ring)
     230            else:
     231                _n = sa2si(base_ring(element),_ring)
     232                _p = p_NSet(_n, _ring)
     233               
     234        # and longs
     235        elif PY_TYPE_CHECK(element, long):
     236            if isinstance(base_ring, FiniteField_prime_modn):
     237                element = element % self.base_ring().characteristic()
     238                _p = p_ISet(int(element),_ring)
     239            else:
     240                _n = sa2si(base_ring(element),_ring)
     241                _p = p_NSet(_n, _ring)
     242
     243        else:
     244            raise NotImplementedError("not able to constructor "+repr(element) +
     245                                      " of type "+ repr(type(element))) #### ??????
     246
     247        return new_NCP(self,_p)
     248
     249
    120250       
    121251    cpdef _coerce_map_from_(self, S):
    122252       """
     
    124254           - the integer ring
    125255           - other localizations away from fewer primes
    126256       """
    127        if S is ZZ:
     257
     258       if self.base_ring().has_coerce_map_from(S):
    128259           return True
    129260       
     261       
     262       
     263    def __hash__(self):
     264       """
     265       Return a hash for this noncommutative ring, that is, a hash of the string
     266       representation of this polynomial ring.
     267
     268       EXAMPLES::
     269           sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     270           sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex')
     271           sage: hash(P)      # somewhat random output
     272           ...
     273
     274       TESTS::
     275
     276       Check conversion from self::
     277           sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     278           sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex')
     279           sage: from sage.matrix.constructor  import Matrix
     280           sage: c = Matrix(3)
     281           sage: c[0,1] = -1
     282           sage: c[0,2] = 1
     283           sage: c[1,2] = 1
     284
     285           sage: from sage.rings.polynomial.plural import NCPolynomialRing_plural
     286           sage: R.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex')
     287           sage: hash(R) == hash(P)
     288           True
     289       """
     290       return hash(self.__repr__())
     291
     292##     def _richcmp_(self, right, int op):
     293##         return (<Parent>self)._richcmp_helper(right, op)
     294     
     295       
     296## ##     cdef int _cmp_c_impl(left, Parent right) except -2:
     297## ##          return 0#cmp(type(left),type(right))
     298
     299##     def __cmp__(self, x):
     300##         """
     301##         EXAMPLES::
     302
     303##         """
     304##         if PY_TYPE_CHECK(x, NCPolynomialRing_plural):
     305##             print "huhu"
     306##             return 0
     307##         print "hihi"
     308##         return 0# return cmp(type(self), type(x))
     309   
     310##     def __richcmp__(self, right, int op):
     311##         print "CMP"
     312##         print "base", self.base_ring()
     313##         print "gens"
     314##         for elt in  self.gens():
     315##             print elt
     316##         print self.term_order()
     317##         for elt in  self.relations():
     318##             print elt
     319##         print  PY_TYPE_CHECK(right, NCPolynomialRing_plural), "huhu"
     320
     321##         print "base", right.base_ring()
     322##         print "gens"
     323##         for elt in  right.gens():
     324##             print elt
     325##         print right.term_order()
     326##         for elt in  right.relations():
     327##             print elt
     328##         return False
     329   
     330##         print cmp(self.base_ring(), right.base_ring())
     331##         print cmp( map(str, self.gens()),  map(str, right.gens()))
     332##         print cmp( map(str, self.relations()),  map(str, right.relations()))
     333
     334##         return False
     335##         if PY_TYPE_CHECK(right, NCPolynomialRing_plural):
     336##             return cmp( (self.base_ring(), map(str, self.gens()),
     337##                          self.term_order(), map(str, self.relations())),
     338##                         (right.base_ring(), map(str, right.gens()),
     339##                          right.term_order(), map(str, right.relations()))
     340##                         )
     341##         else:
     342##             return cmp(type(self),type(right))
     343
     344##     #       return False
     345 
     346 #   def __richcmp__(left, right, int op):
     347 #       return (<ParentWithGens>left)._richcmp(right, op)   
     348 #       return False
     349##         r"""
     350##         Multivariate polynomial rings are said to be equal if:
     351       
     352##         - their base rings match,
     353##         - their generator names match,
     354##         - their term orderings match, and
     355##         - their relations match.
     356
     357
     358##         EXAMPLES::
     359##            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     360##            sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex')
     361
     362##            sage: P == P
     363##            True
     364##            sage: Q = copy(P)
     365##            sage: Q == P
     366##            True
     367                     
     368##            sage: from sage.matrix.constructor  import Matrix
     369##            sage: c = Matrix(3)
     370##            sage: c[0,1] = -1
     371##            sage: c[0,2] = 1
     372##            sage: c[1,2] = 1
     373##            sage: from sage.rings.polynomial.plural import NCPolynomialRing_plural
     374##            sage: R.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex')
     375##            sage: R == P
     376##            True
     377           
     378##            sage: c[0,1] = -2
     379##            sage: R.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex')
     380##            sage: P == R
     381##            False
     382##         """
     383###         return (<Parent>left)._richcmp(right, op)
     384          ##return (<ParentWithGens>left)._richcmp(right, op) 
     385##     cdef int _cmp_c_impl(left, Parent right) except -2:
     386##         r"""
     387##         Multivariate polynomial rings are said to be equal if:
     388       
     389##         - their base rings match,
     390##         - their generator names match,
     391##         - their term orderings match, and
     392##         - their relations match.
     393
     394
     395##         EXAMPLES::
     396##            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     397##            sage: P = A.g_algebra(relations={y*x:-x*y}, order = 'lex')
     398
     399##            sage: P == P
     400##            True
     401##            sage: Q = copy(P)
     402##            sage: Q == P
     403##            True
     404                     
     405##            sage: from sage.matrix.constructor  import Matrix
     406##            sage: c = Matrix(3)
     407##            sage: c[0,1] = -1
     408##            sage: c[0,2] = 1
     409##            sage: c[1,2] = 1
     410##            sage: from sage.rings.polynomial.plural import NCPolynomialRing_plural
     411##            sage: R.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex')
     412##            sage: R == P
     413##            True
     414           
     415##            sage: c[0,1] = -2
     416##            sage: R.<x,y,z> = NCPolynomialRing_plural(QQ, 3, c = c, d = Matrix(3), order='lex')
     417##            sage: P == R
     418##            False
     419##         """       
     420##         print "huhu", PY_TYPE_CHECK(right, NCPolynomialRing_plural)
     421##         if PY_TYPE_CHECK(right, NCPolynomialRing_plural):
     422##             return True
     423##             return cmp( (left.base_ring(), map(str, left.gens()),
     424##                          left.term_order(), map(str, left.relations())),
     425##                         (right.base_ring(), map(str, right.gens()),
     426##                          right.term_order(), map(str, right.relations()))
     427##                         )
     428##         else:
     429##             return cmp(type(left),type(right))
     430       
     431
     432   
    130433    def __pow__(self, n, _):
    131434        """
    132435        Return the free module of rank `n` over this ring.
     
    279582        W = self._list_to_ring(L)
    280583        return new_NRing(W, self.base_ring())
    281584
     585
     586    ### The following methods are handy for implementing Groebner
     587    ### basis algorithms. They do only superficial type/sanity checks
     588    ### and should be called carefully.
     589
     590    def monomial_quotient(self, NCPolynomial_plural f, NCPolynomial_plural g, coeff=False):
     591        r"""
     592        Return ``f/g``, where both ``f`` and`` ``g`` are treated as
     593        monomials.
     594
     595        Coefficients are ignored by default.
     596
     597        INPUT:
     598
     599        - ``f`` - monomial
     600        - ``g`` - monomial
     601        - ``coeff`` - divide coefficients as well (default: ``False``)
     602
     603        EXAMPLES::
     604
     605            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     606            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     607            sage: P.inject_variables()
     608            Defining x, y, z
     609
     610            sage: P.monomial_quotient(3/2*x*y,x,coeff=True)
     611            3/2*y
     612
     613        Note, that `\ZZ` behaves different if ``coeff=True``::
     614
     615            sage: P.monomial_quotient(2*x,3*x)
     616            1
     617
     618            sage: A.<x,y,z> = FreeAlgebra(ZZ, 3)
     619            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     620            sage: P.inject_variables()
     621            Defining x, y, z
     622           
     623            sage: P.monomial_quotient(2*x,3*x,coeff=True)
     624            Traceback (most recent call last):
     625            ...
     626            ArithmeticError: Cannot divide these coefficients.
     627
     628        TESTS::
     629            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     630            sage: R = A.g_algebra(relations={y*x:-x*y},  order='lex')
     631            sage: R.inject_variables()
     632            Defining x, y, z
     633       
     634            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     635            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     636            sage: P.inject_variables()
     637            Defining x, y, z
     638           
     639            sage: P.monomial_quotient(x*y,x)
     640            y
     641
     642##             sage: P.monomial_quotient(x*y,R.gen())
     643##             y
     644
     645            sage: P.monomial_quotient(P(0),P(1))
     646            0
     647
     648            sage: P.monomial_quotient(P(1),P(0))
     649            Traceback (most recent call last):
     650            ...
     651            ZeroDivisionError
     652
     653            sage: P.monomial_quotient(P(3/2),P(2/3), coeff=True)
     654            9/4
     655
     656            sage: P.monomial_quotient(x,y) # Note the wrong result
     657            x*y^1048575*z^1048575 # 64-bit
     658            x*y^65535*z^65535 # 32-bit 
     659
     660            sage: P.monomial_quotient(x,P(1))
     661            x
     662
     663        .. warning::
     664
     665           Assumes that the head term of f is a multiple of the head
     666           term of g and return the multiplicant m. If this rule is
     667           violated, funny things may happen.
     668        """
     669        cdef poly *res
     670        cdef ring *r = self._ring
     671        cdef number *n, *denom
     672       
     673        if not <ParentWithBase>self is f._parent:
     674            f = self._coerce_c(f)
     675        if not <ParentWithBase>self is g._parent:
     676            g = self._coerce_c(g)
     677
     678        if(r != currRing): rChangeCurrRing(r)
     679
     680        if not f._poly:
     681            return self._zero_element
     682        if not g._poly:
     683            raise ZeroDivisionError
     684
     685        res = pDivide(f._poly,g._poly)
     686        if coeff:
     687            if r.ringtype == 0 or r.cf.nDivBy(p_GetCoeff(f._poly, r), p_GetCoeff(g._poly, r)):
     688                n = r.cf.nDiv( p_GetCoeff(f._poly, r) , p_GetCoeff(g._poly, r))
     689                p_SetCoeff0(res, n, r)
     690            else:
     691                raise ArithmeticError("Cannot divide these coefficients.")
     692        else:
     693            p_SetCoeff0(res, n_Init(1, r), r)
     694        return new_NCP(self, res)
     695   
     696    def monomial_divides(self, NCPolynomial_plural a, NCPolynomial_plural b):
     697        """
     698        Return ``False`` if a does not divide b and ``True``
     699        otherwise.
     700
     701        Coefficients are ignored.
     702       
     703        INPUT:
     704
     705        - ``a`` -- monomial
     706
     707        - ``b`` -- monomial
     708
     709        EXAMPLES::
     710
     711            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     712            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     713            sage: P.inject_variables()
     714            Defining x, y, z
     715
     716            sage: P.monomial_divides(x*y*z, x^3*y^2*z^4)
     717            True
     718            sage: P.monomial_divides(x^3*y^2*z^4, x*y*z)
     719            False
     720
     721        TESTS::
     722            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     723            sage: Q = A.g_algebra(relations={y*x:-x*y},  order='lex')
     724            sage: Q.inject_variables()
     725            Defining x, y, z
     726           
     727            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     728            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     729            sage: P.inject_variables()
     730            Defining x, y, z
     731           
     732            sage: P.monomial_divides(P(1), P(0))
     733            True
     734            sage: P.monomial_divides(P(1), x)
     735            True
     736        """
     737        cdef poly *_a
     738        cdef poly *_b
     739        cdef ring *_r
     740        if a._parent is not b._parent:
     741            b = (<NCPolynomialRing_plural>a._parent)._coerce_c(b)
     742
     743        _a = a._poly
     744        _b = b._poly
     745        _r = (<NCPolynomialRing_plural>a._parent)._ring
     746
     747        if _a == NULL:
     748            raise ZeroDivisionError
     749        if _b == NULL:
     750            return True
     751       
     752        if not p_DivisibleBy(_a, _b, _r):
     753            return False
     754        else:
     755            return True
     756
     757
     758    def monomial_lcm(self, NCPolynomial_plural f, NCPolynomial_plural g):
     759        """
     760        LCM for monomials. Coefficients are ignored.
     761       
     762        INPUT:
     763
     764        - ``f`` - monomial
     765       
     766        - ``g`` - monomial
     767
     768        EXAMPLES::
     769
     770            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     771            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     772            sage: P.inject_variables()
     773            Defining x, y, z
     774           
     775            sage: P.monomial_lcm(3/2*x*y,x)
     776            x*y
     777
     778        TESTS::
     779
     780            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     781            sage: R = A.g_algebra(relations={y*x:-x*y},  order='lex')
     782            sage: R.inject_variables()
     783            Defining x, y, z
     784           
     785            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     786            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     787            sage: P.inject_variables()
     788            Defining x, y, z
     789           
     790##             sage: P.monomial_lcm(x*y,R.gen())
     791##             x*y
     792
     793            sage: P.monomial_lcm(P(3/2),P(2/3))
     794            1
     795
     796            sage: P.monomial_lcm(x,P(1))
     797            x
     798        """
     799        cdef poly *m = p_ISet(1,self._ring)
     800       
     801        if not <ParentWithBase>self is f._parent:
     802            f = self._coerce_c(f)
     803        if not <ParentWithBase>self is g._parent:
     804            g = self._coerce_c(g)
     805
     806        if f._poly == NULL:
     807            if g._poly == NULL:
     808                return self._zero_element
     809            else:
     810                raise ArithmeticError, "Cannot compute LCM of zero and nonzero element."
     811        if g._poly == NULL:
     812            raise ArithmeticError, "Cannot compute LCM of zero and nonzero element."
     813
     814        if(self._ring != currRing): rChangeCurrRing(self._ring)
     815       
     816        pLcm(f._poly, g._poly, m)
     817        p_Setm(m, self._ring)
     818        return new_NCP(self,m)
     819       
     820    def monomial_reduce(self, NCPolynomial_plural f, G):
     821        """
     822        Try to find a ``g`` in ``G`` where ``g.lm()`` divides
     823        ``f``. If found ``(flt,g)`` is returned, ``(0,0)`` otherwise,
     824        where ``flt`` is ``f/g.lm()``.
     825
     826        It is assumed that ``G`` is iterable and contains *only*
     827        elements in this polynomial ring.
     828
     829        Coefficients are ignored.
     830       
     831        INPUT:
     832
     833        - ``f`` - monomial
     834        - ``G`` - list/set of mpolynomials
     835           
     836        EXAMPLES::
     837
     838            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     839            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     840            sage: P.inject_variables()
     841            Defining x, y, z
     842
     843            sage: f = x*y^2
     844            sage: G = [ 3/2*x^3 + y^2 + 1/2, 1/4*x*y + 2/7, 1/2  ]
     845            sage: P.monomial_reduce(f,G)
     846            (y, 1/4*x*y + 2/7)
     847
     848        TESTS::
     849            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     850            sage: Q = A.g_algebra(relations={y*x:-x*y},  order='lex')
     851            sage: Q.inject_variables()
     852            Defining x, y, z
     853           
     854            sage: A.<x,y,z> = FreeAlgebra(QQ, 3)
     855            sage: P = A.g_algebra(relations={y*x:-x*y},  order='lex')
     856            sage: P.inject_variables()
     857            Defining x, y, z
     858            sage: f = x*y^2
     859            sage: G = [ 3/2*x^3 + y^2 + 1/2, 1/4*x*y + 2/7, 1/2  ]
     860
     861            sage: P.monomial_reduce(P(0),G)
     862            (0, 0)
     863
     864            sage: P.monomial_reduce(f,[P(0)])
     865            (0, 0)
     866        """
     867        cdef poly *m = f._poly
     868        cdef ring *r = self._ring
     869        cdef poly *flt
     870
     871        if not m:
     872            return f,f
     873       
     874        for g in G:
     875            if PY_TYPE_CHECK(g, NCPolynomial_plural) \
     876                   and (<NCPolynomial_plural>g) \
     877                   and p_LmDivisibleBy((<NCPolynomial_plural>g)._poly, m, r):
     878                flt = pDivide(f._poly, (<NCPolynomial_plural>g)._poly)
     879                #p_SetCoeff(flt, n_Div( p_GetCoeff(f._poly, r) , p_GetCoeff((<NCPolynomial_plural>g)._poly, r), r), r)
     880                p_SetCoeff(flt, n_Init(1, r), r)
     881                return new_NCP(self,flt), g
     882        return self._zero_element,self._zero_element
     883
     884    def monomial_pairwise_prime(self, NCPolynomial_plural g, NCPolynomial_plural h):
     885        """
     886        Return ``True`` if ``h`` and ``g`` are pairwise prime. Both
     887        are treated as monomials.
     888
     889        Coefficients are ignored.
     890
     891        INPUT:
     892
     893        - ``h`` - monomial
     894        - ``g`` - monomial
     895
     896        EXAMPLES::
     897
     898            sage: P.<x,y,z> = QQ[]
     899            sage: P.monomial_pairwise_prime(x^2*z^3, y^4)
     900            True
     901
     902            sage: P.monomial_pairwise_prime(1/2*x^3*y^2, 3/4*y^3)
     903            False
     904
     905        TESTS::
     906
     907            sage: Q.<x,y,z> = QQ[]
     908            sage: P.<x,y,z> = QQ[]
     909            sage: P.monomial_pairwise_prime(x^2*z^3, Q('y^4'))
     910            True
     911
     912            sage: P.monomial_pairwise_prime(1/2*x^3*y^2, Q(0))
     913            True
     914
     915            sage: P.monomial_pairwise_prime(P(1/2),x)
     916            False
     917        """
     918        cdef int i
     919        cdef ring *r
     920        cdef poly *p, *q
     921
     922        if h._parent is not g._parent:
     923            g = (<MPolynomialRing_libsingular>h._parent)._coerce_c(g)
     924
     925        r = (<MPolynomialRing_libsingular>h._parent)._ring
     926        p = g._poly
     927        q = h._poly
     928
     929        if p == NULL:
     930            if q == NULL:
     931                return False #GCD(0,0) = 0
     932            else:
     933                return True #GCD(x,0) = 1
     934
     935        elif q == NULL:
     936            return True # GCD(0,x) = 1
     937
     938        elif p_IsConstant(p,r) or p_IsConstant(q,r): # assuming a base field
     939            return False
     940
     941        for i from 1 <= i <= r.N:
     942            if p_GetExp(p,i,r) and p_GetExp(q,i,r):
     943                return False
     944        return True
     945
     946    def monomial_all_divisors(self, NCPolynomial_plural t):
     947        """
     948        Return a list of all monomials that divide ``t``.
     949
     950        Coefficients are ignored.
     951         
     952        INPUT:
     953
     954        - ``t`` - a monomial
     955 
     956        OUTPUT:
     957            a list of monomials
     958
     959
     960        EXAMPLES::
     961
     962            sage: P.<x,y,z> = QQ[]
     963            sage: P.monomial_all_divisors(x^2*z^3)
     964            [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]
     965           
     966        ALGORITHM: addwithcarry idea by Toon Segers
     967        """
     968
     969        M = list()
     970
     971        cdef ring *_ring = self._ring
     972        cdef poly *maxvector = t._poly
     973        cdef poly *tempvector = p_ISet(1, _ring)
     974         
     975        pos = 1
     976         
     977        while not p_ExpVectorEqual(tempvector, maxvector, _ring):
     978          tempvector = addwithcarry(tempvector, maxvector, pos, _ring)
     979          M.append(new_NCP(self, p_Copy(tempvector,_ring)))
     980        return M
     981
     982
     983
    282984cdef class NCPolynomial_plural(RingElement):
    283985    def __init__(self, NCPolynomialRing_plural parent):
    284986        self._poly = NULL
     
    9201622    I = H.ideal([H.gen(i) *H.gen(i) for i in alt_vars]).twostd()
    9211623    return H.quotient(I)
    9221624
     1625cdef poly *addwithcarry(poly *tempvector, poly *maxvector, int pos, ring *_ring):
     1626    if p_GetExp(tempvector, pos, _ring) < p_GetExp(maxvector, pos, _ring):
     1627      p_SetExp(tempvector, pos, p_GetExp(tempvector, pos, _ring)+1, _ring)
     1628    else:
     1629      p_SetExp(tempvector, pos, 0, _ring)
     1630      tempvector = addwithcarry(tempvector, maxvector, pos + 1, _ring)
     1631    p_Setm(tempvector, _ring)
     1632    return tempvector