Ticket #4539: extplural-more2.patch

File extplural-more2.patch, 24.5 KB (added by AlexanderDreyer, 4 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