Ticket #7644: trac_7644_reversion_lagrange_4.patch

File trac_7644_reversion_lagrange_4.patch, 18.6 KB (added by Niles Johnson, 12 years ago)

apply only this patch; choice for default precision improved

  • doc/en/reference/power_series.rst

    # HG changeset patch
    # User Niles Johnson <nilesj@gmail.com>
    # Date 1278638817 14400
    # Node ID 4e3101381089da2be01407ef78c49c47feda6632
    # Parent  ce0a8069b40f05035bdfc5b6831c91212e2f4919
    [mq]: trac_7644_reversion_lagrange_3.patch
    
    diff -r ce0a8069b40f -r 4e3101381089 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/laurent_series_ring
    13    sage/rings/laurent_series_ring_element
    14  No newline at end of file
     14   sage/rings/laurent_series_ring_element
  • sage/rings/power_series_poly.pyx

    diff -r ce0a8069b40f -r 4e3101381089 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, then the precision
     707        of the reversion defaults to the default precision of
     708        ``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, and self must have finite precision (i.e. this cannot be done
     712        for polynomials).  Under the current implementation, the leading
     713        coefficient must be a unit in the base ring, and the base ring must
     714        have characteristic zero.
    665715
    666         EXAMPLES:
     716        ALGORITHM::
     717
     718        We first attempt to pass the computation to pari; if this fails, we
     719        use Lagrange inversion.  Using ``sage: set_verbose(1)`` will print
     720        a message if passing to pari fails.
     721
     722        EXAMPLES::
     723       
    667724            sage: R.<x> = PowerSeriesRing(QQ)
    668             sage: f = 2*x + 3*x**2 - x**4 + O(x**5)
     725            sage: f = 2*x + 3*x^2 - x^4 + O(x^5)
    669726            sage: g = f.reversion()
    670727            sage: g
    671728            1/2*x - 3/8*x^2 + 9/16*x^3 - 131/128*x^4 + O(x^5)
     
    674731            sage: g(f)
    675732            x + O(x^5)
    676733
    677             sage: f += 1
     734            sage: A.<t> = PowerSeriesRing(ZZ)
     735            sage: a = t - t^2 - 2*t^4 + t^5 + O(t^6)
     736            sage: b = a.reversion(); b
     737            t + t^2 + 2*t^3 + 7*t^4 + 25*t^5 + O(t^6)
     738            sage: a(b)
     739            t + O(t^6)
     740            sage: b(a)
     741            t + O(t^6)
     742
     743        If the leading coefficient is not a unit, we pass to its fraction
     744        field if possible::
     745                 
     746            sage: A.<t> = PowerSeriesRing(ZZ)
     747            sage: a = 2*t - 4*t^2 + t^4 - t^5 + O(t^6)
     748            sage: a.reversion()
     749            1/2*t + 1/2*t^2 + t^3 + 79/32*t^4 + 437/64*t^5 + O(t^6)
     750       
     751            sage: B.<b> = PolynomialRing(ZZ)
     752            sage: A.<t> = PowerSeriesRing(B)
     753            sage: f = 2*b*t + b*t^2 + 3*b^2*t^3 + O(t^4)
     754            sage: g = f.reversion(); g
     755            1/(2*b)*t - 1/(8*b^2)*t^2 + ((-3*b + 1)/(16*b^3))*t^3 + O(t^4)
     756            sage: g.base_ring()
     757            Fraction Field of Univariate Polynomial Ring in b over Integer Ring
     758            sage: f(g)
     759            t + O(t^4)
     760            sage: g(f)
     761            t + O(t^4)
     762
     763        In positive characteristic, attempt first to lift to characteristic
     764        zero and perform the reversion there::
     765       
     766            sage: A8.<t> = PowerSeriesRing(Zmod(8))
     767            sage: a = t - 15*t^2 - 2*t^4 + t^5 + O(t^6)
     768            sage: b = a.reversion(); b
     769            t + 7*t^2 + 2*t^3 + 5*t^4 + t^5 + O(t^6)
     770            sage: a(b)
     771            t + O(t^6)
     772            sage: b(a)
     773            t + O(t^6)
     774
     775
     776        The optional argument ``precision`` sets the precision of the output::
     777       
     778            sage: R.<x> = PowerSeriesRing(QQ)
     779            sage: f = 2*x + 3*x^2 - 7*x^3 + x^4 + O(x^5)
     780            sage: g = f.reversion(precision=3); g
     781            1/2*x - 3/8*x^2 + O(x^3)
     782            sage: f(g)
     783            x + O(x^3)
     784            sage: g(f)
     785            x + O(x^3)
     786
     787        If the input series has infinite precision, the precision of the
     788        output is automatically set to the default precision of the parent
     789        ring::
     790
     791            sage: R.<x> = PowerSeriesRing(QQ, default_prec=20)
     792            sage: (x - x^2).reversion() # get some Catalan numbers
     793            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)
     794            sage: (x - x^2).reversion(precision=3)
     795            x + x^2 + O(x^3)
     796
     797
     798        TESTS::
     799
     800            sage: R.<x> = PowerSeriesRing(QQ)
     801            sage: f = 1 + 2*x + 3*x**2 - x**4 + O(x**5)
    678802            sage: f.reversion()
    679803            Traceback (most recent call last):
    680804            ...
    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
     805            ValueError: Series must have valuation one for reversion.
     806
     807
     808
    686809        """
    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"
    691810        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())
     811            raise ValueError("Series must have valuation one for reversion.")
     812
     813        f = self
     814
     815        if f.prec() is infinity and precision is None:
     816            precision = f.parent().default_prec()
     817        if precision:
     818            f = f.add_bigoh(precision)
     819
     820        out_prec = f.prec()
     821
     822        if not f[1].is_unit():
     823            # if leading coefficient is not a unit, attempt passing
     824            # to fraction field
     825            try:
     826                f = f.change_ring(f.base_ring().fraction_field())
     827            except TypeError:
     828                raise TypeError("Leading coefficient must be a unit, or base ring must have a fraction field.")
     829
     830        # set output parent after possibly passing to fraction field,
     831        # but before possibly lifting to characteristic zero
     832        out_parent = f.parent()
     833
     834        # first, try reversion with pari; this is faster than Lagrange inversion
     835        try:
     836            f2 = f._pari_()
     837            g = f2.serreverse()
     838            return PowerSeries_poly(f.parent(),g.Vecrev(),out_prec)
     839        except (TypeError,ValueError,AttributeError):
     840            # if pari fails, continue with Lagrange inversion
     841            from sage.misc.all import verbose
     842            verbose("passing to pari failed; trying Lagrange inversion")
     843           
     844       
     845        if f.parent().characteristic() > 0:
     846            # over a ring of positive characteristic, attempt lifting to
     847            # characteristic zero ring
     848            verbose("parent ring has positive characteristic; attempting lift to characteristic zero")
     849            base_lift = f.base_ring().lift().codomain()
     850            verbose("characteristic zero base is "+str(base_lift))
     851            f_lift = f.change_ring(base_lift)
     852            verbose("f_lift is "+str(f_lift))
     853            rev_lift = f_lift.reversion()
     854            return rev_lift.change_ring(f.base_ring())
     855       
     856        t = f.parent().gen()
     857
     858        h = t/f
     859        k = 1
     860        g = 0
     861        for i in range(1, out_prec):
     862            k *= h
     863            g += k.padded_list(i)[i - 1]/i*t**i
     864        g = g.add_bigoh(out_prec)
     865        return PowerSeries_poly(out_parent, g, out_prec, check=False)
     866
    696867
    697868def make_powerseries_poly_v0(parent,  f, prec, is_gen):
    698869    """
     
    703874    instead make a new function and make sure that both kinds of
    704875    objects correctly unpickle as the new type.
    705876
    706     EXAMPLES:
     877    EXAMPLES::
     878   
    707879        sage: R.<t> = QQ[[]]
    708880        sage: sage.rings.power_series_poly.make_powerseries_poly_v0(R, t, infinity, True)
    709881        t