Ticket #7644: trac_7644_reversion_lagrange_6.patch

File trac_7644_reversion_lagrange_6.patch, 18.8 KB (added by Niles Johnson, 12 years ago)

apply only this patch

  • doc/en/reference/power_series.rst

    # HG changeset patch
    # User Niles Johnson <nilesj@gmail.com>
    # Date 1278638817 14400
    # Node ID 7866809677caacfa0585d2b83e28da42d6c09580
    # Parent  8597fede2103fb473914d49eab97a4886016141c
    Ticket 7644 univariate power series reversion using Lagrange inversion
    
    diff -r 8597fede2103 -r 7866809677ca doc/en/reference/power_series.rst
    a b  
    88   
    99   sage/rings/power_series_ring
    1010   sage/rings/power_series_ring_element
     11   sage/rings/power_series_poly
    1112
    1213   sage/rings/multi_power_series_ring
    1314   sage/rings/multi_power_series_ring_element
  • sage/rings/power_series_poly.pyx

    diff -r 8597fede2103 -r 7866809677ca sage/rings/power_series_poly.pyx
    a b  
     1"""
     2Power Series Methods
     3
     4The class ``PowerSeries_poly`` provides additional methods for univariate power series.
     5"""
     6
     7
    18include "../ext/stdsage.pxi"
    29
    310from power_series_ring_element cimport PowerSeries
     
    1219
    1320    def __init__(self, parent, f=0, prec=infinity, int check=1, is_gen=0):
    1421        """
    15         EXAMPLES:
     22        EXAMPLES::
     23       
    1624            sage: R, q = PowerSeriesRing(CC, 'q').objgen()
    1725            sage: R
    1826            Power Series Ring in q over Complex Field with 53 bits of precision
     
    5159        """
    5260        Return a hash of self.
    5361       
    54         EXAMPLES:
     62        EXAMPLES::
     63       
    5564            sage: R.<t> = ZZ[[]]
    5665            sage: t.__hash__()
    5766            760233507         # 32-bit
     
    6675        """
    6776        Used for pickling.
    6877
    69         EXAMPLES:
     78        EXAMPLES::
     79       
    7080            sage: A.<z> = RR[[]]
    7181            sage: f = z - z^3 + O(z^10)
    7282            sage: f == loads(dumps(f)) # uses __reduce__
     
    7989       """
    8090       Used for comparing power series.
    8191
    82        EXAMPLES:
     92       EXAMPLES::
     93       
    8394           sage: R.<t> = ZZ[[]]
    8495           sage: f = 1 + t + t^7 - 5*t^10
    8596           sage: g = 1 + t + t^7 - 5*t^10 + O(t^15)
     
    94105       
    95106    def polynomial(self):
    96107        """
    97         EXAMPLE:
     108        Return the underlying polynomial of self.
     109
     110        EXAMPLE::
     111       
    98112            sage: R.<t> = GF(7)[[]]
    99113            sage: f = 3 - t^3 + O(t^5)
    100114            sage: f.polynomial()
     
    106120        """
    107121        Return the valuation of self.
    108122
    109         EXAMPLES:
     123        EXAMPLES::
     124       
    110125            sage: R.<t> = QQ[[]]
    111126            sage: (5 - t^8 + O(t^11)).valuation()
    112127            0
     
    121136       
    122137    def degree(self):
    123138        """
    124         Return the degree of the polynomial associated to self. That
     139        Return the degree of the underlying polynomial of self. That
    125140        is, if self is of the form f(x) + O(x^n), we return the degree
    126141        of f(x). Note that if f(x) is 0, we return -1, just as with
    127142        polynomials.
    128143
    129         EXAMPLES:
     144        EXAMPLES::
     145       
    130146            sage: R.<t> = ZZ[[]]
    131147            sage: (5 + t^3 + O(t^4)).degree()
    132148            3
     
    141157        """
    142158        Return True if self is nonzero, and False otherwise.
    143159
    144         EXAMPLES:
     160        EXAMPLES::
     161       
    145162            sage: R.<t> = GF(11)[[]]
    146163            sage: (1 + t + O(t^18)).__nonzero__()
    147164            True
     
    154171
    155172    def __call__(self, *xs):
    156173        """
    157         EXAMPLE:
     174        EXAMPLES::
     175       
    158176            sage: R.<t> = GF(7)[[]]
    159177            sage: f = 3 - t^3 + O(t^5)
    160178            sage: f(1)
     
    193211
    194212        ** DO NOT USE THIS ** unless you know what you're doing.
    195213
    196         EXAMPLES:
     214        EXAMPLES::
     215       
    197216            sage: R.<t> = GF(7)[[]]
    198217            sage: f = 3 + 6*t^3 + O(t^5)
    199218            sage: f._unsafe_mutate(0, 5)
     
    202221            sage: f._unsafe_mutate(2, 1) ; f
    203222            5 + t^2 + 6*t^3 + O(t^5)
    204223           
    205         Mutating can even bump up the precision.
     224        - Mutating can even bump up the precision::
     225       
    206226            sage: f._unsafe_mutate(6, 1) ; f
    207227            5 + t^2 + 6*t^3 + t^6 + O(t^7)
    208228            sage: f._unsafe_mutate(0, 0) ; f
     
    230250        Returns 0 for negative coefficients. Raises an IndexError if
    231251        try to access beyond known coefficients.
    232252
    233         EXAMPLES:
     253        EXAMPLES::
     254       
    234255            sage: R.<t> = QQ[[]]
    235256            sage: f = 3/2 - 17/5*t^3 + O(t^5)
    236257            sage: f[3]
     
    313334        """
    314335        Return an iterator over the coefficients of this power series.
    315336
    316         EXAMPLES:
     337        EXAMPLES::
     338       
    317339            sage: R.<t> = QQ[[]]
    318340            sage: f = t + 17/5*t^3 + 2*t^4 + O(t^5)
    319341            sage: for a in f: print a,
     
    325347        """
    326348        Return the negative of this power series.
    327349
    328         EXAMPLES:
     350        EXAMPLES::
     351       
    329352            sage: R.<t> = QQ[[]]
    330353            sage: f = t + 17/5*t^3 + 2*t^4 + O(t^5)
    331354            sage: -f
     
    336359       
    337360    cpdef ModuleElement _add_(self, ModuleElement right_m):
    338361        """
    339         EXAMPLES:
     362        EXAMPLES::
     363       
    340364            sage: R.<x> = PowerSeriesRing(ZZ)
    341365            sage: f = x^4 + O(x^5); f
    342366            x^4 + O(x^5)
     
    351375                                         
    352376    cpdef ModuleElement _iadd_(self, ModuleElement right_m):
    353377        """
    354         EXAMPLES:
     378        EXAMPLES::
     379       
    355380            sage: R.<x> = PowerSeriesRing(ZZ)
    356381            sage: f = x^4
    357382            sage: f += x; f
     
    380405        """
    381406        Return the difference of two power series.
    382407
    383         EXAMPLES:
     408        EXAMPLES::
     409       
    384410            sage: k.<w> = ZZ[]
    385411            sage: R.<t> = k[[]]
    386412            sage: w*t^2 -w*t +13 - (w*t^2 + w*t)
     
    394420        """
    395421        Return the product of two power series.
    396422
    397         EXAMPLES:
     423        EXAMPLES::
     424       
    398425            sage: k.<w> = ZZ[[]]
    399426            sage: (1+17*w+15*w^3+O(w^5))*(19*w^10+O(w^12))
    400427            19*w^10 + 323*w^11 + O(w^12)       
     
    409436        """
    410437        Set self to self * right_r, and return this result.
    411438
    412         EXAMPLES:
     439        EXAMPLES::
     440       
    413441            sage: k.<w> = ZZ[[]]
    414442            sage: f = (1+17*w+15*w^3+O(w^5))
    415443            sage: f *= (19*w^10+O(w^12))
     
    431459        """
    432460        Multiply self on the right by a scalar.
    433461
    434         EXAMPLES:
     462        EXAMPLES::
     463       
    435464            sage: R.<t> = GF(7)[[]]
    436465            sage: f = t + 3*t^4 + O(t^11)
    437466            sage: f * GF(7)(3)
     
    443472        """
    444473        Multiply self on the left by a scalar.
    445474
    446         EXAMPLES:
     475        EXAMPLES::
     476       
    447477            sage: R.<t> = GF(11)[[]]
    448478            sage: f = 1 + 3*t^4 + O(t^120)
    449479            sage: 2 * f
     
    455485        """
    456486        Set self to self left-multiplied by a scalar.
    457487
    458         EXAMPLES:
     488        EXAMPLES::
     489       
    459490            sage: R.<t> = GF(13)[[]]
    460491            sage: f = 3 + 7*t^3 + O(t^4)
    461492            sage: f._ilmul_(2)
     
    468499   
    469500    def __floordiv__(self, denom):
    470501        """
    471         EXAMPLES:
     502        EXAMPLES::
     503       
    472504            sage: R.<t> = ZZ[[]] ; f = t**10-1 ; g = 1+t+t^7 ; h = f.add_bigoh(20)
    473505            sage: f // g
    474506            -1 + t - t^2 + t^3 - t^4 + t^5 - t^6 + 2*t^7 - 3*t^8 + 4*t^9 - 4*t^10 + 5*t^11 - 6*t^12 + 7*t^13 - 9*t^14 + 12*t^15 - 16*t^16 + 20*t^17 - 25*t^18 + 31*t^19 + O(t^20)
     
    497529        """
    498530        Shift self to the left by n, i.e. multiply by x^n.
    499531
    500         EXAMPLES:
     532        EXAMPLES::
     533       
    501534            sage: R.<t> = QQ[[]]
    502535            sage: f = 1 + t + t^4
    503536            sage: f << 1
     
    513546        Shift self to the right by n, i.e. multiply by x^-n and
    514547        remove any terms of negative exponent.
    515548
    516         EXAMPLES:
     549        EXAMPLES::
     550       
    517551            sage: R.<t> = GF(2)[[]]
    518552            sage: f = t + t^4 + O(t^7)
    519553            sage: f >> 1
     
    528562
    529563    def truncate(self, prec=infinity):
    530564        """
    531         The polynomial obtained from power series by truncation.
     565        The polynomial obtained from power series by truncation at
     566        precision ``prec``.
    532567
    533         EXAMPLES:
     568        EXAMPLES::
     569       
    534570            sage: R.<I> = GF(2)[[]]
    535571            sage: f = 1/(1+I+O(I^8)); f
    536572            1 + I + I^2 + I^3 + I^4 + I^5 + I^6 + I^7 + O(I^8)
     
    544580           
    545581    cdef _inplace_truncate(self, long prec):
    546582        """
    547         Truncate self to precision prec in place.
     583        Truncate self to precision ``prec`` in place.
    548584
    549         NOTE: This is very unsafe, since power series are supposed to
    550         be immutable in Sage. Use at your own risk!
     585        NOTE::
     586
     587            This is very unsafe, since power series are supposed to
     588            be immutable in Sage. Use at your own risk!
    551589        """
    552590        self.__f = self.__f._inplace_truncate(prec)
    553591        self.prec = prec
     
    555593
    556594    def truncate_powerseries(self, long prec):
    557595        r"""
    558         Returns the power series of degree $ < n$ which is equivalent to self
    559         modulo $x^n$.
     596        Given input ``prec`` = $n$, returns the power series of degree
     597        $< n$ which is equivalent to self modulo $x^n$.
    560598
    561         EXAMPLES:
     599        EXAMPLES::
     600       
    562601            sage: R.<I> = GF(2)[[]]
    563602            sage: f = 1/(1+I+O(I^8)); f
    564603            1 + I + I^2 + I^3 + I^4 + I^5 + I^6 + I^7 + O(I^8)
     
    574613        the list of coefficients of the underlying polynomial, so in
    575614        particular, need not have length equal to self.prec().
    576615
    577         EXAMPLES:
     616        EXAMPLES::
     617       
    578618            sage: R.<t> = ZZ[[]]
    579619            sage: f = 1 - 5*t^3 + t^5 + O(t^7)
    580620            sage: f.list()
     
    588628        dict for the underlying polynomial, so need not have keys
    589629        corresponding to every number smaller than self.prec().
    590630
    591         EXAMPLES:
     631        EXAMPLES::
     632       
    592633            sage: R.<t> = ZZ[[]]
    593634            sage: f = 1 + t^10 + O(t^12)
    594635            sage: f.dict()
     
    607648        Otherwise, we call _derivative(var) on each coefficient of
    608649        the series.
    609650
    610         SEE ALSO:
     651        SEE ALSO::
     652       
    611653            self.derivative()
    612654       
    613         EXAMPLES:
     655        EXAMPLES::
     656       
    614657            sage: R.<t> = PowerSeriesRing(QQ, sparse=True)
    615658            sage: f = 2 + 3*t^2 + t^100000 + O(t^10000000); f
    616659            2 + 3*t^2 + t^100000 + O(t^10000000)
     
    642685        """
    643686        The integral of this power series with 0 constant term.
    644687
    645         EXAMPLES:
     688        EXAMPLES::
     689       
    646690            sage: k.<w> = QQ[[]]
    647691            sage: (1+17*w+15*w^3+O(w^5)).integral()
    648692            w + 17/2*w^2 + 15/4*w^4 + O(w^6)       
     
    654698        return PowerSeries_poly(self._parent, self.__f.integral(),
    655699                                         self.prec()+1, check=False)
    656700
    657     def reversion(self):
     701    def reversion(self, precision=None):
    658702        """
    659         Return the reversion of f, i.e., the series g such that
    660         g(f(x)) = x.
     703        Return the reversion of f, i.e., the series g such that g(f(x)) =
     704        x.  Given an optional argument ``precision``, return the reversion
     705        with given precision (note that the reversion can have precision at
     706        most ``f.prec()``).  If ``f`` has infinite precision, and the argument
     707        ``precision`` is not given, then the precision of the reversion
     708        defaults to the default precision of ``f.parent()``.
    661709
    662         Note that this is only possible if self.valuation() is exactly
    663         1, and must have finite precision (i.e. this cannot be done
    664         for polynomials).
     710        Note that this is only possible if the valuation of self is exactly
     711        1.
     712       
     713        ALGORITHM:
    665714
    666         EXAMPLES:
     715        We first attempt to pass the computation to pari; if this fails, we
     716        use Lagrange inversion.  Using ``sage: set_verbose(1)`` will print
     717        a message if passing to pari fails.
     718
     719        If the base ring has positive characteristic, then we attempt to
     720        lift to a characteristic zero ring and perform the reversion there.
     721        If this fails, an error is raised.
     722
     723        EXAMPLES::
     724       
    667725            sage: R.<x> = PowerSeriesRing(QQ)
    668             sage: f = 2*x + 3*x**2 - x**4 + O(x**5)
     726            sage: f = 2*x + 3*x^2 - x^4 + O(x^5)
    669727            sage: g = f.reversion()
    670728            sage: g
    671729            1/2*x - 3/8*x^2 + 9/16*x^3 - 131/128*x^4 + O(x^5)
     
    674732            sage: g(f)
    675733            x + O(x^5)
    676734
    677             sage: f += 1
     735            sage: A.<t> = PowerSeriesRing(ZZ)
     736            sage: a = t - t^2 - 2*t^4 + t^5 + O(t^6)
     737            sage: b = a.reversion(); b
     738            t + t^2 + 2*t^3 + 7*t^4 + 25*t^5 + O(t^6)
     739            sage: a(b)
     740            t + O(t^6)
     741            sage: b(a)
     742            t + O(t^6)
     743
     744            sage: B.<b,c> = PolynomialRing(ZZ)
     745            sage: A.<t> = PowerSeriesRing(B)
     746            sage: f = t + b*t^2 + c*t^3 + O(t^4)
     747            sage: g = f.reversion(); g
     748            t - b*t^2 + (2*b^2 - c)*t^3 + O(t^4)
     749            sage: f(g)
     750            t + O(t^4)
     751            sage: g(f)
     752            t + O(t^4)
     753
     754
     755        If the leading coefficient is not a unit, we pass to its fraction
     756        field if possible::
     757                 
     758            sage: A.<t> = PowerSeriesRing(ZZ)
     759            sage: a = 2*t - 4*t^2 + t^4 - t^5 + O(t^6)
     760            sage: a.reversion()
     761            1/2*t + 1/2*t^2 + t^3 + 79/32*t^4 + 437/64*t^5 + O(t^6)
     762       
     763            sage: B.<b> = PolynomialRing(ZZ)
     764            sage: A.<t> = PowerSeriesRing(B)
     765            sage: f = 2*b*t + b*t^2 + 3*b^2*t^3 + O(t^4)
     766            sage: g = f.reversion(); g
     767            1/(2*b)*t - 1/(8*b^2)*t^2 + ((-3*b + 1)/(16*b^3))*t^3 + O(t^4)
     768            sage: g.base_ring()
     769            Fraction Field of Univariate Polynomial Ring in b over Integer Ring
     770            sage: f(g)
     771            t + O(t^4)
     772            sage: g(f)
     773            t + O(t^4)
     774
     775        We can handle some base rings of positive characteristic::
     776       
     777            sage: A8.<t> = PowerSeriesRing(Zmod(8))
     778            sage: a = t - 15*t^2 - 2*t^4 + t^5 + O(t^6)
     779            sage: b = a.reversion(); b
     780            t + 7*t^2 + 2*t^3 + 5*t^4 + t^5 + O(t^6)
     781            sage: a(b)
     782            t + O(t^6)
     783            sage: b(a)
     784            t + O(t^6)
     785
     786        The optional argument ``precision`` sets the precision of the output::
     787       
     788            sage: R.<x> = PowerSeriesRing(QQ)
     789            sage: f = 2*x + 3*x^2 - 7*x^3 + x^4 + O(x^5)
     790            sage: g = f.reversion(precision=3); g
     791            1/2*x - 3/8*x^2 + O(x^3)
     792            sage: f(g)
     793            x + O(x^3)
     794            sage: g(f)
     795            x + O(x^3)
     796
     797        If the input series has infinite precision, the precision of the
     798        output is automatically set to the default precision of the parent
     799        ring::
     800
     801            sage: R.<x> = PowerSeriesRing(QQ, default_prec=20)
     802            sage: (x - x^2).reversion() # get some Catalan numbers
     803            x + x^2 + 2*x^3 + 5*x^4 + 14*x^5 + 42*x^6 + 132*x^7 + 429*x^8 + 1430*x^9 + 4862*x^10 + 16796*x^11 + 58786*x^12 + 208012*x^13 + 742900*x^14 + 2674440*x^15 + 9694845*x^16 + 35357670*x^17 + 129644790*x^18 + 477638700*x^19 + O(x^20)
     804            sage: (x - x^2).reversion(precision=3)
     805            x + x^2 + O(x^3)
     806
     807
     808        TESTS::
     809
     810            sage: R.<x> = PowerSeriesRing(QQ)
     811            sage: f = 1 + 2*x + 3*x^2 - x^4 + O(x^5)
    678812            sage: f.reversion()
    679813            Traceback (most recent call last):
    680814            ...
    681             ValueError: series must have valuation one for reversion
    682             sage: x.reversion()
    683             Traceback (most recent call last):
    684             ...
    685             ValueError: series must have finite precision for reversion
     815            ValueError: Series must have valuation one for reversion.
     816
     817
     818
    686819        """
    687         if not isinstance(self.parent().base_ring(), rational_field.RationalField):
    688             raise NotImplementedError
    689         if self.prec() is infinity:
    690             raise ValueError, "series must have finite precision for reversion"
    691820        if self.valuation() != 1:
    692             raise ValueError, "series must have valuation one for reversion"
    693         f = self._pari_()
    694         g = f.serreverse()
    695         return PowerSeries_poly(self.parent(),g.Vecrev(),self.prec())
     821            raise ValueError("Series must have valuation one for reversion.")
     822
     823        f = self
     824
     825        if f.prec() is infinity and precision is None:
     826            precision = f.parent().default_prec()
     827        if precision:
     828            f = f.add_bigoh(precision)
     829
     830        out_prec = f.prec()
     831
     832        if not f[1].is_unit():
     833            # if leading coefficient is not a unit, attempt passing
     834            # to fraction field
     835            try:
     836                f = f.change_ring(f.base_ring().fraction_field())
     837            except TypeError:
     838                raise TypeError("Leading coefficient must be a unit, or base ring must have a fraction field.")
     839
     840        # set output parent after possibly passing to fraction field,
     841        # but before possibly lifting to characteristic zero
     842        out_parent = f.parent()
     843
     844        # first, try reversion with pari; this is faster than Lagrange inversion
     845        try:
     846            f2 = f._pari_()
     847            g = f2.serreverse()
     848            return PowerSeries_poly(f.parent(),g.Vecrev(),out_prec)
     849        except (TypeError,ValueError,AttributeError):
     850            # if pari fails, continue with Lagrange inversion
     851            from sage.misc.all import verbose
     852            verbose("passing to pari failed; trying Lagrange inversion")
     853           
     854       
     855        if f.parent().characteristic() > 0:
     856            # over a ring of positive characteristic, attempt lifting to
     857            # characteristic zero ring
     858            verbose("parent ring has positive characteristic; attempting lift to characteristic zero")
     859            base_lift = f.base_ring().lift().codomain()
     860            verbose("characteristic zero base is "+str(base_lift))
     861            f_lift = f.change_ring(base_lift)
     862            verbose("f_lift is "+str(f_lift))
     863            rev_lift = f_lift.reversion()
     864            return rev_lift.change_ring(f.base_ring())
     865       
     866        t = f.parent().gen()
     867
     868        h = t/f
     869        k = 1
     870        g = 0
     871        for i in range(1, out_prec):
     872            k *= h
     873            g += k.padded_list(i)[i - 1]/i*t**i
     874        g = g.add_bigoh(out_prec)
     875        return PowerSeries_poly(out_parent, g, out_prec, check=False)
     876
    696877
    697878def make_powerseries_poly_v0(parent,  f, prec, is_gen):
    698879    """
     
    703884    instead make a new function and make sure that both kinds of
    704885    objects correctly unpickle as the new type.
    705886
    706     EXAMPLES:
     887    EXAMPLES::
     888   
    707889        sage: R.<t> = QQ[[]]
    708890        sage: sage.rings.power_series_poly.make_powerseries_poly_v0(R, t, infinity, True)
    709891        t