Ticket #7644: trac_7644_reversion_lagrange_6.patch
File trac_7644_reversion_lagrange_6.patch, 18.8 KB (added by , 12 years ago) 


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 8 8 9 9 sage/rings/power_series_ring 10 10 sage/rings/power_series_ring_element 11 sage/rings/power_series_poly 11 12 12 13 sage/rings/multi_power_series_ring 13 14 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 """ 2 Power Series Methods 3 4 The class ``PowerSeries_poly`` provides additional methods for univariate power series. 5 """ 6 7 1 8 include "../ext/stdsage.pxi" 2 9 3 10 from power_series_ring_element cimport PowerSeries … … 12 19 13 20 def __init__(self, parent, f=0, prec=infinity, int check=1, is_gen=0): 14 21 """ 15 EXAMPLES: 22 EXAMPLES:: 23 16 24 sage: R, q = PowerSeriesRing(CC, 'q').objgen() 17 25 sage: R 18 26 Power Series Ring in q over Complex Field with 53 bits of precision … … 51 59 """ 52 60 Return a hash of self. 53 61 54 EXAMPLES: 62 EXAMPLES:: 63 55 64 sage: R.<t> = ZZ[[]] 56 65 sage: t.__hash__() 57 66 760233507 # 32bit … … 66 75 """ 67 76 Used for pickling. 68 77 69 EXAMPLES: 78 EXAMPLES:: 79 70 80 sage: A.<z> = RR[[]] 71 81 sage: f = z  z^3 + O(z^10) 72 82 sage: f == loads(dumps(f)) # uses __reduce__ … … 79 89 """ 80 90 Used for comparing power series. 81 91 82 EXAMPLES: 92 EXAMPLES:: 93 83 94 sage: R.<t> = ZZ[[]] 84 95 sage: f = 1 + t + t^7  5*t^10 85 96 sage: g = 1 + t + t^7  5*t^10 + O(t^15) … … 94 105 95 106 def polynomial(self): 96 107 """ 97 EXAMPLE: 108 Return the underlying polynomial of self. 109 110 EXAMPLE:: 111 98 112 sage: R.<t> = GF(7)[[]] 99 113 sage: f = 3  t^3 + O(t^5) 100 114 sage: f.polynomial() … … 106 120 """ 107 121 Return the valuation of self. 108 122 109 EXAMPLES: 123 EXAMPLES:: 124 110 125 sage: R.<t> = QQ[[]] 111 126 sage: (5  t^8 + O(t^11)).valuation() 112 127 0 … … 121 136 122 137 def degree(self): 123 138 """ 124 Return the degree of the polynomial associated toself. That139 Return the degree of the underlying polynomial of self. That 125 140 is, if self is of the form f(x) + O(x^n), we return the degree 126 141 of f(x). Note that if f(x) is 0, we return 1, just as with 127 142 polynomials. 128 143 129 EXAMPLES: 144 EXAMPLES:: 145 130 146 sage: R.<t> = ZZ[[]] 131 147 sage: (5 + t^3 + O(t^4)).degree() 132 148 3 … … 141 157 """ 142 158 Return True if self is nonzero, and False otherwise. 143 159 144 EXAMPLES: 160 EXAMPLES:: 161 145 162 sage: R.<t> = GF(11)[[]] 146 163 sage: (1 + t + O(t^18)).__nonzero__() 147 164 True … … 154 171 155 172 def __call__(self, *xs): 156 173 """ 157 EXAMPLE: 174 EXAMPLES:: 175 158 176 sage: R.<t> = GF(7)[[]] 159 177 sage: f = 3  t^3 + O(t^5) 160 178 sage: f(1) … … 193 211 194 212 ** DO NOT USE THIS ** unless you know what you're doing. 195 213 196 EXAMPLES: 214 EXAMPLES:: 215 197 216 sage: R.<t> = GF(7)[[]] 198 217 sage: f = 3 + 6*t^3 + O(t^5) 199 218 sage: f._unsafe_mutate(0, 5) … … 202 221 sage: f._unsafe_mutate(2, 1) ; f 203 222 5 + t^2 + 6*t^3 + O(t^5) 204 223 205 Mutating can even bump up the precision. 224  Mutating can even bump up the precision:: 225 206 226 sage: f._unsafe_mutate(6, 1) ; f 207 227 5 + t^2 + 6*t^3 + t^6 + O(t^7) 208 228 sage: f._unsafe_mutate(0, 0) ; f … … 230 250 Returns 0 for negative coefficients. Raises an IndexError if 231 251 try to access beyond known coefficients. 232 252 233 EXAMPLES: 253 EXAMPLES:: 254 234 255 sage: R.<t> = QQ[[]] 235 256 sage: f = 3/2  17/5*t^3 + O(t^5) 236 257 sage: f[3] … … 313 334 """ 314 335 Return an iterator over the coefficients of this power series. 315 336 316 EXAMPLES: 337 EXAMPLES:: 338 317 339 sage: R.<t> = QQ[[]] 318 340 sage: f = t + 17/5*t^3 + 2*t^4 + O(t^5) 319 341 sage: for a in f: print a, … … 325 347 """ 326 348 Return the negative of this power series. 327 349 328 EXAMPLES: 350 EXAMPLES:: 351 329 352 sage: R.<t> = QQ[[]] 330 353 sage: f = t + 17/5*t^3 + 2*t^4 + O(t^5) 331 354 sage: f … … 336 359 337 360 cpdef ModuleElement _add_(self, ModuleElement right_m): 338 361 """ 339 EXAMPLES: 362 EXAMPLES:: 363 340 364 sage: R.<x> = PowerSeriesRing(ZZ) 341 365 sage: f = x^4 + O(x^5); f 342 366 x^4 + O(x^5) … … 351 375 352 376 cpdef ModuleElement _iadd_(self, ModuleElement right_m): 353 377 """ 354 EXAMPLES: 378 EXAMPLES:: 379 355 380 sage: R.<x> = PowerSeriesRing(ZZ) 356 381 sage: f = x^4 357 382 sage: f += x; f … … 380 405 """ 381 406 Return the difference of two power series. 382 407 383 EXAMPLES: 408 EXAMPLES:: 409 384 410 sage: k.<w> = ZZ[] 385 411 sage: R.<t> = k[[]] 386 412 sage: w*t^2 w*t +13  (w*t^2 + w*t) … … 394 420 """ 395 421 Return the product of two power series. 396 422 397 EXAMPLES: 423 EXAMPLES:: 424 398 425 sage: k.<w> = ZZ[[]] 399 426 sage: (1+17*w+15*w^3+O(w^5))*(19*w^10+O(w^12)) 400 427 19*w^10 + 323*w^11 + O(w^12) … … 409 436 """ 410 437 Set self to self * right_r, and return this result. 411 438 412 EXAMPLES: 439 EXAMPLES:: 440 413 441 sage: k.<w> = ZZ[[]] 414 442 sage: f = (1+17*w+15*w^3+O(w^5)) 415 443 sage: f *= (19*w^10+O(w^12)) … … 431 459 """ 432 460 Multiply self on the right by a scalar. 433 461 434 EXAMPLES: 462 EXAMPLES:: 463 435 464 sage: R.<t> = GF(7)[[]] 436 465 sage: f = t + 3*t^4 + O(t^11) 437 466 sage: f * GF(7)(3) … … 443 472 """ 444 473 Multiply self on the left by a scalar. 445 474 446 EXAMPLES: 475 EXAMPLES:: 476 447 477 sage: R.<t> = GF(11)[[]] 448 478 sage: f = 1 + 3*t^4 + O(t^120) 449 479 sage: 2 * f … … 455 485 """ 456 486 Set self to self leftmultiplied by a scalar. 457 487 458 EXAMPLES: 488 EXAMPLES:: 489 459 490 sage: R.<t> = GF(13)[[]] 460 491 sage: f = 3 + 7*t^3 + O(t^4) 461 492 sage: f._ilmul_(2) … … 468 499 469 500 def __floordiv__(self, denom): 470 501 """ 471 EXAMPLES: 502 EXAMPLES:: 503 472 504 sage: R.<t> = ZZ[[]] ; f = t**101 ; g = 1+t+t^7 ; h = f.add_bigoh(20) 473 505 sage: f // g 474 506 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) … … 497 529 """ 498 530 Shift self to the left by n, i.e. multiply by x^n. 499 531 500 EXAMPLES: 532 EXAMPLES:: 533 501 534 sage: R.<t> = QQ[[]] 502 535 sage: f = 1 + t + t^4 503 536 sage: f << 1 … … 513 546 Shift self to the right by n, i.e. multiply by x^n and 514 547 remove any terms of negative exponent. 515 548 516 EXAMPLES: 549 EXAMPLES:: 550 517 551 sage: R.<t> = GF(2)[[]] 518 552 sage: f = t + t^4 + O(t^7) 519 553 sage: f >> 1 … … 528 562 529 563 def truncate(self, prec=infinity): 530 564 """ 531 The polynomial obtained from power series by truncation. 565 The polynomial obtained from power series by truncation at 566 precision ``prec``. 532 567 533 EXAMPLES: 568 EXAMPLES:: 569 534 570 sage: R.<I> = GF(2)[[]] 535 571 sage: f = 1/(1+I+O(I^8)); f 536 572 1 + I + I^2 + I^3 + I^4 + I^5 + I^6 + I^7 + O(I^8) … … 544 580 545 581 cdef _inplace_truncate(self, long prec): 546 582 """ 547 Truncate self to precision precin place.583 Truncate self to precision ``prec`` in place. 548 584 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! 551 589 """ 552 590 self.__f = self.__f._inplace_truncate(prec) 553 591 self.prec = prec … … 555 593 556 594 def truncate_powerseries(self, long prec): 557 595 r""" 558 Returns the power series of degree $ < n$ which is equivalent to self559 modulo $x^n$.596 Given input ``prec`` = $n$, returns the power series of degree 597 $< n$ which is equivalent to self modulo $x^n$. 560 598 561 EXAMPLES: 599 EXAMPLES:: 600 562 601 sage: R.<I> = GF(2)[[]] 563 602 sage: f = 1/(1+I+O(I^8)); f 564 603 1 + I + I^2 + I^3 + I^4 + I^5 + I^6 + I^7 + O(I^8) … … 574 613 the list of coefficients of the underlying polynomial, so in 575 614 particular, need not have length equal to self.prec(). 576 615 577 EXAMPLES: 616 EXAMPLES:: 617 578 618 sage: R.<t> = ZZ[[]] 579 619 sage: f = 1  5*t^3 + t^5 + O(t^7) 580 620 sage: f.list() … … 588 628 dict for the underlying polynomial, so need not have keys 589 629 corresponding to every number smaller than self.prec(). 590 630 591 EXAMPLES: 631 EXAMPLES:: 632 592 633 sage: R.<t> = ZZ[[]] 593 634 sage: f = 1 + t^10 + O(t^12) 594 635 sage: f.dict() … … 607 648 Otherwise, we call _derivative(var) on each coefficient of 608 649 the series. 609 650 610 SEE ALSO: 651 SEE ALSO:: 652 611 653 self.derivative() 612 654 613 EXAMPLES: 655 EXAMPLES:: 656 614 657 sage: R.<t> = PowerSeriesRing(QQ, sparse=True) 615 658 sage: f = 2 + 3*t^2 + t^100000 + O(t^10000000); f 616 659 2 + 3*t^2 + t^100000 + O(t^10000000) … … 642 685 """ 643 686 The integral of this power series with 0 constant term. 644 687 645 EXAMPLES: 688 EXAMPLES:: 689 646 690 sage: k.<w> = QQ[[]] 647 691 sage: (1+17*w+15*w^3+O(w^5)).integral() 648 692 w + 17/2*w^2 + 15/4*w^4 + O(w^6) … … 654 698 return PowerSeries_poly(self._parent, self.__f.integral(), 655 699 self.prec()+1, check=False) 656 700 657 def reversion(self ):701 def reversion(self, precision=None): 658 702 """ 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()``. 661 709 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: 665 714 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 667 725 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) 669 727 sage: g = f.reversion() 670 728 sage: g 671 729 1/2*x  3/8*x^2 + 9/16*x^3  131/128*x^4 + O(x^5) … … 674 732 sage: g(f) 675 733 x + O(x^5) 676 734 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) 678 812 sage: f.reversion() 679 813 Traceback (most recent call last): 680 814 ... 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 686 819 """ 687 if not isinstance(self.parent().base_ring(), rational_field.RationalField):688 raise NotImplementedError689 if self.prec() is infinity:690 raise ValueError, "series must have finite precision for reversion"691 820 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 696 877 697 878 def make_powerseries_poly_v0(parent, f, prec, is_gen): 698 879 """ … … 703 884 instead make a new function and make sure that both kinds of 704 885 objects correctly unpickle as the new type. 705 886 706 EXAMPLES: 887 EXAMPLES:: 888 707 889 sage: R.<t> = QQ[[]] 708 890 sage: sage.rings.power_series_poly.make_powerseries_poly_v0(R, t, infinity, True) 709 891 t