Ticket #9944: trac9944_polynomial_speedup.patch

File trac9944_polynomial_speedup.patch, 16.5 KB (added by SimonKing, 8 years ago)

Speedup of polynomial arithmetic by improved base ring conversion

  • sage/categories/algebras.py

    # HG changeset patch
    # User Simon King <simon.king@uni-jena.de>
    # Date 1302703896 -7200
    # Node ID e48ba7940e6a9f0abca53621c48c8b7eccf30687
    # Parent  3f897822f54ddce8138c5a7cc9b8b392f4f991b2
    #9944: Speedup conversion from base ring. Call Parent.__init__ in ring initialization.
    
    diff --git a/sage/categories/algebras.py b/sage/categories/algebras.py
    a b  
    9898
    9999        def __init_extra__(self):
    100100            """
    101             Declares the canonical coercion from ``self.base_ring()`` to ``self``.
     101            Declares the canonical coercion from ``self.base_ring()`` to ``self``,
     102            if there has been none before.
    102103
    103104            EXAMPLES::
    104105
     
    112113                sage: A(1)          # indirect doctest
    113114                B[word: ]
    114115            """
    115             # TODO: find and implement a good syntax!
    116             self.register_coercion(SetMorphism(function = self.from_base_ring, parent = Hom(self.base_ring(), self, Rings())))
    117116            # Should be a morphism of Algebras(self.base_ring()), but e.g. QQ is not in Algebras(QQ)
     117            mor = SetMorphism(function = self.from_base_ring, parent = Hom(self.base_ring(), self, Rings()))
     118            try:
     119                self.register_coercion(mor)
     120            except AssertionError:
     121                pass
    118122
    119123    class ElementMethods:
    120124        # TODO: move the content of AlgebraElement here
  • sage/rings/polynomial/multi_polynomial_libsingular.pyx

    diff --git a/sage/rings/polynomial/multi_polynomial_libsingular.pyx b/sage/rings/polynomial/multi_polynomial_libsingular.pyx
    a b  
    3333- Martin Albrecht (2009-06): refactored the code to allow better
    3434  re-use
    3535
     36- Simon King (2011-03): Use a faster way of conversion from the base
     37  ring.
     38
    3639TODO:
    3740
    3841- implement Real, Complex coefficient rings via libSINGULAR
     
    293296            sage: P.<x,y,z> = PolynomialRing(Integers(2^32),order='lex')
    294297            sage: P(2^32-1)
    295298            4294967295
     299
     300        TEST:
     301
     302        Make sure that a faster conversion map from the base ring is used;
     303        see trac ticket #9944::
     304
     305            sage: R.<x,y> = PolynomialRing(ZZ)
     306            sage: R.convert_map_from(R.base_ring())
     307            Polynomial base injection morphism:
     308              From: Integer Ring
     309              To:   Multivariate Polynomial Ring in x, y over Integer Ring
     310
    296311        """
    297312        MPolynomialRing_generic.__init__(self, base_ring, n, names, order)
    298313        self._has_singular = True
     
    301316        self._ring = singular_ring_new(base_ring, n, self._names, order)
    302317        self._one_element = <MPolynomial_libsingular>new_MP(self,p_ISet(1, self._ring))
    303318        self._zero_element = <MPolynomial_libsingular>new_MP(self,NULL)
     319        # This polynomial ring should belong to Algebras(base_ring).
     320        # Algebras(...).parent_class, which was called from MPolynomialRing_generic.__init__,
     321        # tries to provide a conversion from the base ring, if it does not exist.
     322        # This is for algebras that only do the generic stuff in their initialisation.
     323        # But here, we want to use PolynomialBaseringInjection. Hence, we need to
     324        # wipe the memory and construct the conversion from scratch.
     325        from sage.rings.polynomial.polynomial_element import PolynomialBaseringInjection
     326        base_inject = PolynomialBaseringInjection(base_ring, self)
     327        self.register_conversion(base_inject)
    304328
    305329    def __dealloc__(self):
    306330        r"""
  • sage/rings/polynomial/multi_polynomial_ring.py

    diff --git a/sage/rings/polynomial/multi_polynomial_ring.py b/sage/rings/polynomial/multi_polynomial_ring.py
    a b  
    150150        self._gens = tuple(self._gens)
    151151        self._zero_tuple = tuple(v)
    152152        self._has_singular = can_convert_to_singular(self)
     153        # This polynomial ring should belong to Algebras(base_ring).
     154        # Algebras(...).parent_class, which was called from MPolynomialRing_generic.__init__,
     155        # tries to provide a conversion from the base ring, if it does not exist.
     156        # This is for algebras that only do the generic stuff in their initialisation.
     157        # But here, we want to use PolynomialBaseringInjection. Hence, we need to
     158        # wipe the memory and construct the conversion from scratch.
     159        if n:
     160            from sage.rings.polynomial.polynomial_element import PolynomialBaseringInjection
     161            base_inject = PolynomialBaseringInjection(base_ring, self)
     162            self.register_conversion(base_inject)
    153163
    154164    def _monomial_order_function(self):
    155165        return self.__monomial_order_function
  • sage/rings/polynomial/polynomial_element.pyx

    diff --git a/sage/rings/polynomial/polynomial_element.pyx b/sage/rings/polynomial/polynomial_element.pyx
    a b  
    33
    44AUTHORS:
    55
    6 -  William Stein: first version
     6-  William Stein: first version.
    77
    88-  Martin Albrecht: Added singular coercion.
    99
    10 -  Robert Bradshaw: Move Polynomial_generic_dense to Cython
     10-  Robert Bradshaw: Move Polynomial_generic_dense to Cython.
     11
     12-  Simon King: Use a faster way of conversion from the base ring.
    1113
    1214
    1315TESTS::
     
    60206022    This class is used for conversion from a ring to a polynomial
    60216023    over that ring.
    60226024   
    6023     It calls the _new_constant_poly method on the generator,
    6024     which should be optimized for a particular polynomial type.
     6025    If the polynomial ring has a One and if the elements provide an
     6026    ``_rmul_`` method, and ``_rmul_`` does *not* return None (which is
     6027    the case for `p`-adics) then conversion is obtained by multiplying
     6028    the base ring element with the One by means of ``_rmul_``.
     6029
     6030    Otherwise It calls the _new_constant_poly method on the generator,
     6031    which should be optimized for a particular polynomial type, but
     6032    often isn't.
     6033
    60256034    Technically, it should be a method of the polynomial ring, but
    60266035    few polynomial rings are cython classes.
     6036
     6037    NOTE:
     6038
     6039    We use ``_rmul_`` and not ``_lmul_`` since for many polynomial rings
     6040    ``_lmul_`` simply calls ``_rmul_``.
     6041
     6042    EXAMPLES:
     6043
     6044    We demonstrate that different polynomial ring classes use
     6045    polynomial base injection maps::
     6046
     6047        sage: R.<x> = Qp(3)[]
     6048        sage: R.convert_map_from(R.base_ring())
     6049        Polynomial base injection morphism:
     6050          From: 3-adic Field with capped relative precision 20
     6051          To:   Univariate Polynomial Ring in x over 3-adic Field with capped relative precision 20
     6052        sage: R.<x,y> = Qp(3)[]
     6053        sage: R.convert_map_from(R.base_ring())
     6054        Polynomial base injection morphism:
     6055          From: 3-adic Field with capped relative precision 20
     6056          To:   Multivariate Polynomial Ring in x, y over 3-adic Field with capped relative precision 20
     6057        sage: R.<x,y> = QQ[]
     6058        sage: R.convert_map_from(R.base_ring())
     6059        Polynomial base injection morphism:
     6060          From: Rational Field
     6061          To:   Multivariate Polynomial Ring in x, y over Rational Field
     6062        sage: R.<x> = QQ[]
     6063        sage: R.convert_map_from(R.base_ring())
     6064        Polynomial base injection morphism:
     6065          From: Rational Field
     6066          To:   Univariate Polynomial Ring in x over Rational Field
     6067
     6068    By trac ticket #9944, there are now only very few exceptions::
     6069
     6070        sage: PolynomialRing(QQ,names=[]).convert_map_from(QQ)
     6071        Call morphism:
     6072          From: Rational Field
     6073          To:   Multivariate Polynomial Ring in no variables over Rational Field
     6074
    60276075    """
    60286076   
    6029     cdef Polynomial _an_element
     6077    cdef RingElement _an_element
     6078    cdef object _one_rmul_
    60306079   
    60316080    def __init__(self, domain, codomain):
    60326081        """
    6033         TESTS:
     6082        TESTS::
     6083
    60346084            sage: from sage.rings.polynomial.polynomial_element import PolynomialBaseringInjection
    60356085            sage: PolynomialBaseringInjection(QQ, QQ['x'])
    60366086            Polynomial base injection morphism:
     
    60406090            Traceback (most recent call last):
    60416091            ...
    60426092            AssertionError: domain must be basering
     6093
     6094        ::
     6095
     6096            sage: R.<t> = Qp(2)[]
     6097            sage: f = R.convert_map_from(R.base_ring())    # indirect doctest
     6098            sage: f(Qp(2).one()*3)
     6099            (1 + 2 + O(2^20))
     6100            sage: (Qp(2).one()*3)*t
     6101            (1 + 2 + O(2^20))*t
     6102
    60436103        """
    60446104        assert domain is codomain.base_ring(), "domain must be basering"
    60456105        Morphism.__init__(self, domain, codomain)
    60466106        self._an_element = codomain.gen()
    60476107        self._repr_type_str = "Polynomial base injection"
     6108        if domain is codomain: # some rings are base rings of themselves!
     6109            return
     6110        try:
     6111            one = codomain._element_constructor_(domain.one_element())
     6112        except (AttributeError, NotImplementedError, TypeError):
     6113            # perhaps it uses the old model?
     6114            try:
     6115                one = codomain._coerce_c(domain.one_element())
     6116            except (AttributeError, NotImplementedError, TypeError):
     6117                return
     6118        try:
     6119            one_rmul_ = one._rmul_
     6120        except AttributeError:
     6121            return
     6122        # For the p-adic fields, _lmul_ and _rmul_ return None!!!
     6123        # To work around, we need to test its sanity before we try
     6124        # to use it.
     6125        try:
     6126            if one_rmul_(domain.one_element()) is None:
     6127                return
     6128            self._one_rmul_ = one_rmul_
     6129        except TypeError:
     6130            pass
    60486131   
    60496132    cpdef Element _call_(self, x):
    60506133        """
     
    60546137            Polynomial base injection morphism:
    60556138              From: Integer Ring
    60566139              To:   Univariate Polynomial Ring in x over Integer Ring
    6057             sage: m(2) # implicit doctest
     6140            sage: m(2) # indirect doctest
    60586141            2
    60596142            sage: parent(m(2))
    60606143            Univariate Polynomial Ring in x over Integer Ring
    60616144        """
     6145        if self._one_rmul_ is not None:
     6146            return self._one_rmul_(x)
    60626147        return self._an_element._new_constant_poly(<Element>x)
    60636148   
    60646149    cpdef Element _call_with_args(self, x, args=(), kwds={}):
     
    60666151        TESTS:
    60676152            sage: from sage.rings.polynomial.polynomial_element import PolynomialBaseringInjection
    60686153            sage: m = PolynomialBaseringInjection(Qp(5), Qp(5)['x'])
    6069             sage: m(1 + O(5^11), absprec = 5)
     6154            sage: m(1 + O(5^11), absprec = 5)   # indirect doctest
    60706155            (1 + O(5^11))
    60716156        """
    60726157        return self._codomain._element_constructor(x, *args, **kwds)
  • sage/rings/polynomial/polynomial_element_generic.py

    diff --git a/sage/rings/polynomial/polynomial_element_generic.py b/sage/rings/polynomial/polynomial_element_generic.py
    a b  
    460460        output.__normalize()
    461461        return output
    462462
     463    def _rmul_(self, left):
     464        r"""
     465        EXAMPLES::
     466
     467            sage: R.<x> = PolynomialRing(ZZ, sparse=True)
     468            sage: (x^100000 - x^50000) * (x^100000 + x^50000)
     469            x^200000 - x^100000
     470            sage: 7 * (x^100000 - x^50000)   # indirect doctest
     471            7*x^100000 - 7*x^50000
     472
     473        AUTHOR:
     474
     475        - Simon King (2011-03-31)
     476        """
     477        output = {}
     478
     479        for (index, coeff) in self.__coeffs.iteritems():
     480            output[index] = left * coeff
     481
     482        output = self.parent()(output, check=False)
     483        output.__normalize()
     484        return output
     485
     486    def _lmul_(self, right):
     487        r"""
     488        EXAMPLES::
     489
     490            sage: R.<x> = PolynomialRing(ZZ, sparse=True)
     491            sage: (x^100000 - x^50000) * (x^100000 + x^50000)
     492            x^200000 - x^100000
     493            sage: (x^100000 - x^50000) * 7   # indirect doctest
     494            7*x^100000 - 7*x^50000
     495
     496        AUTHOR:
     497
     498        - Simon King (2011-03-31)
     499        """
     500        output = {}
     501
     502        for (index, coeff) in self.__coeffs.iteritems():
     503            output[index] = coeff * right
     504
     505        output = self.parent()(output, check=False)
     506        output.__normalize()
     507        return output
     508
    463509    def shift(self, n):
    464510        r"""
    465511        Returns this polynomial multiplied by the power `x^n`. If `n` is negative,
  • sage/rings/polynomial/polynomial_ring.py

    diff --git a/sage/rings/polynomial/polynomial_ring.py b/sage/rings/polynomial/polynomial_ring.py
    a b  
    199199            category = categories.IntegralDomains()
    200200        else:
    201201            category = categories.CommutativeRings()
    202         sage.algebras.algebra.Algebra.__init__(self, base_ring, names=name, normalize=True, category=category)
    203202        self.__is_sparse = sparse
    204203        if element_class:
    205204            self._polynomial_class = element_class
     
    209208            else:
    210209                from sage.rings.polynomial import polynomial_element
    211210                self._polynomial_class = polynomial_element.Polynomial_generic_dense
    212         self.__generator = self._polynomial_class(self, [0,1], is_gen=True)
    213211        self.__cyclopoly_cache = {}
    214212        self._has_singular = False
     213        # Algebra.__init__ also calls __init_extra__ of Algebras(...).parent_class, which
     214        # tries to provide a conversion from the base ring, if it does not exist.
     215        # This is for algebras that only do the generic stuff in their initialisation.
     216        # But here, we want to use PolynomialBaseringInjection. Hence, we need to
     217        # wipe the memory and construct the conversion from scratch.
     218        sage.algebras.algebra.Algebra.__init__(self, base_ring, names=name, normalize=True, category=category)
     219        self.__generator = self._polynomial_class(self, [0,1], is_gen=True)
     220        base_inject = PolynomialBaseringInjection(base_ring, self)
     221        self._unset_coercions_used()
    215222        self._populate_coercion_lists_(
    216                 coerce_list = [PolynomialBaseringInjection(base_ring, self)],
    217                 convert_list = [list],
     223                coerce_list = [base_inject],
     224                convert_list = [list, base_inject],
    218225                convert_method_name = '_polynomial_')
    219226 
    220227
  • sage/rings/qqbar.py

    diff --git a/sage/rings/qqbar.py b/sage/rings/qqbar.py
    a b  
    375375    R.<x> = QQbar[]
    376376    v1 = AA(2)
    377377    v2 = QQbar(sqrt(v1))
    378     v3 = QQbar(3)
    379     v4 = sqrt(v3)
    380     v5 = v2*v4
    381     v6 = (1 - v2)*(1 - v4) - 1 - v5
    382     v7 = QQbar(sqrt(v1))
    383     v8 = sqrt(v3)
    384     si1 = v7*v8
    385     cp = AA.common_polynomial(x^2 + ((1 - v7)*(1 + v8) - 1 + si1)*x - si1)
    386     v9 = QQbar.polynomial_root(cp, RIF(-RR(1.7320508075688774), -RR(1.7320508075688772)))
    387     v10 = 1 - v9
    388     v11 = v6 + (v10 - 1)
    389     v12 = -1 - v4 - QQbar.polynomial_root(cp, RIF(-RR(1.7320508075688774), -RR(1.7320508075688772)))
    390     v13 = 1 + v12
    391     v14 = v10*(v6 + v5) - (v6 - v5*v9)
    392     si2 = v5*v9
    393     AA.polynomial_root(AA.common_polynomial(x^4 + (v11 + (v13 - 1))*x^3 + (v14 + (v13*v11 - v11))*x^2 + (v13*(v14 - si2) - (v14 - si2*v12))*x - si2*v12), RIF(RR(0.99999999999999989), RR(1.0000000000000002)))
     378    v3 = sqrt(QQbar(3))
     379    v4 = v2*v3
     380    v5 = (1 - v2)*(1 - v3) - 1 - v4
     381    v6 = QQbar(sqrt(v1))
     382    si1 = v6*v3
     383    cp = AA.common_polynomial(x^2 + ((1 - v6)*(1 + v3) - 1 + si1)*x - si1)
     384    v7 = QQbar.polynomial_root(cp, RIF(-RR(1.7320508075688774), -RR(1.7320508075688772)))
     385    v8 = 1 - v7
     386    v9 = v5 + (v8 - 1)
     387    v10 = -1 - v3 - QQbar.polynomial_root(cp, RIF(-RR(1.7320508075688774), -RR(1.7320508075688772)))
     388    v11 = 1 + v10
     389    v12 = v8*(v5 + v4) - (v5 - v4*v7)
     390    si2 = v4*v7
     391    AA.polynomial_root(AA.common_polynomial(x^4 + (v9 + (v11 - 1))*x^3 + (v12 + (v11*v9 - v9))*x^2 + (v11*(v12 - si2) - (v12 - si2*v10))*x - si2*v10), RIF(RR(0.99999999999999989), RR(1.0000000000000002)))
    394392    sage: one
    395393    1
    396394
  • sage/rings/ring.pyx

    diff --git a/sage/rings/ring.pyx b/sage/rings/ring.pyx
    a b  
    7272    """
    7373    Generic ring class.
    7474    """
     75    # Unfortunately, ParentWithGens inherits from sage.structure.parent_old.Parent.
     76    # Its __init__ method does *not* call Parent.__init__, since this would somehow
     77    # yield an infinite recursion. But when we call it from here, it works.
     78    # This is done in order to ensure that __init_extra__ is called.
     79    def __init__(self, base, names=None, normalize=True, category = None):
     80        ParentWithGens.__init__(self, base, names=names, normalize=normalize)
     81        Parent.__init__(self, category=category)
     82
    7583    def __iter__(self):
    7684        r"""
    7785        Return an iterator through the elements of self. Not implemented in general.