Ticket #10532: trac_10532_invert.patch

File trac_10532_invert.patch, 2.4 KB (added by pernici, 9 years ago)
  • sage/rings/multi_power_series_ring_element.py

    # HG changeset patch
    # User Mario Pernici <mario.pernici@gmail.com>
    # Date 1297188793 -3600
    # Node ID 34aa7eff3b5a026162349dc6658d6a75a4f8fb54
    # Parent  e0d02abc6dca0215ea783fc42404684403e1a3b2
    In this patch __invert__ is implemented using _mul_trunc_generic;
    inversion of a generic multivariate series becomes much faster.
    
    diff -r e0d02abc6dca -r 34aa7eff3b5a sage/rings/multi_power_series_ring_element.py
    a b  
    174174from sage.structure.category_object import CategoryObject
    175175from sage.rings.polynomial.polynomial_element import Polynomial_generic_dense
    176176from sage.rings.power_series_poly import PowerSeries_poly
    177 
     177import sage.misc.misc
    178178
    179179def is_MPowerSeries(f):
    180180    """
     
    563563            1 - a - b + a^2 + 3*a*b + a*c + b^2 + b*c - a^3 - 5*a^2*b
    564564            - 2*a^2*c - 5*a*b^2 - 4*a*b*c - b^3 - 2*b^2*c + O(a, b, c)^4
    565565        """
    566         if self.valuation() == 0:
    567             return self.parent(self._bg_value.__invert__())
    568         else:
     566        if self.valuation() != 0:
    569567            raise NotImplementedError("Multiplicative inverse of multivariate power series currently implemented only if constant coefficient is a unit.")
     568        if self == 1:
     569            return self
     570        s = self._bg_value
     571        prec = s.prec()
     572        if prec is infinity and s.degree() > 0:
     573            prec = s._parent.default_prec()
     574
     575        # Use Newton's method, i.e. start with single term approximation,
     576        # and then iteratively compute $x' = 2x - Ax^2$, where $A$ is the
     577        # series we're trying to invert.
     578       
     579        try:
     580            first_coeff = ~s[0]
     581        except ValueError, ZeroDivisionError:
     582            raise ZeroDivisionError, "leading coefficient must be a unit"
     583
     584        if prec is infinity:
     585            return s._parent(first_coeff, prec=prec)
     586
     587        A = s.truncate()
     588        R = A.parent()     # R is the corresponding polynomial ring
     589        current = R(first_coeff)
     590        for next_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
     591            z = current.square()._mul_trunc_generic(A.truncate(next_prec),next_prec)
     592            current = 2*current - z
     593        return self._parent(s.parent()(current, prec=prec))
    570594
    571595    ## comparisons
    572596    def __cmp__(self, other):