Ticket #9706: trac_9706_chebyshev_patch_new.patch
File trac_9706_chebyshev_patch_new.patch, 14.5 KB (added by , 8 years ago) 


sage/functions/orthogonal_polys.py
# HG changeset patch # User Stefan Reiterer <domors@gmx.net> # Date 1386270617 3600 # Node ID 92488b626f7f9a296ee7e0d83550dfd8eac57c73 # Parent 12e04e4ed442d2a6b6d9af7919f87ee49191776d Trac 9706: Added Chebyshev Polynomials, and Base class. Major corrections and Bugfixes. diff git a/sage/functions/orthogonal_polys.py b/sage/functions/orthogonal_polys.py
a b 287 287 288 288  :wikipedia:`Associated_Legendre_polynomials` 289 289 290  [EffCheby] Wolfram Koepf: Effcient Computation of Chebyshev Polynomials 291 in Computer Algebra 292 Computer Algebra Systems: A Practical Guide. 293 John Wiley, Chichester (1999): 7999. 290 294 291 295 AUTHORS: 292 296 … … 448 452 """ 449 453 raise ValueError("no special values known") 450 454 451 def __call__(self,n,x):452 """453 This overides the call method from SageObject to avoid problems with coercions,454 since the _eval_ method is able to handle for more data types than symbolic functions455 would normally allow.456 457 EXAMPLES::458 459 sage: K.<a> = NumberField(x^3x1)460 sage: chebyshev_T(5, a)461 16*a^2 + a  4462 463 464 """465 try:466 return super(OrthogonalPolynomial,self).__call__(n,x)467 except TypeError:468 return self._eval_(n,x)469 470 471 455 def _eval_(self, *args): 472 456 """ 473 457 The symbolic evaluation is either done with maxima, or with direct … … 487 471 sage: chebyshev_T(n,1) 488 472 (1)^n 489 473 sage: chebyshev_T(7,x) 490 chebyshev_T(7, x)474 64*x^7  112*x^5 + 56*x^3  7*x 491 475 sage: chebyshev_T(3/2,x) 492 476 chebyshev_T(3/2, x) 493 477 """ … … 495 479 if not is_Expression(args[1]) and is_inexact(args[1]): 496 480 try: 497 481 import sage.libs.mpmath.all as mpmath 498 father = parent(args[1]) # convert back to original type 499 return father(self._evalf_(*args)) 482 return self._evalf_(*args) 500 483 except AttributeError: 501 484 pass 502 485 except mpmath.NoConvergence: … … 514 497 pass 515 498 else: 516 499 return None 517 518 if args[0] < 0:519 return None520 521 500 try: 522 501 return self._eval_special_values_(*args) 523 502 except ValueError: 524 503 pass 525 504 505 #if negative indices are not specified 506 #in _eval_special_values only return symbolic 507 #value 508 if args[0] < 0 and args[0] in ZZ: 509 return None 510 526 511 if not is_Expression(args[0]): 527 512 try: 528 513 return self._clenshaw_method_(*args) … … 544 529 545 530 return s.sage() 546 531 532 def __call__(self,n,x): 533 """ 534 This overides the call method from SageObject to avoid problems with coercions, 535 since the _eval_ method is able to handle for more data types than symbolic functions 536 would normally allow. 537 Thus we have the distinction between algebraic objects (if n is an integer), 538 and else as symbolic function. 539 540 EXAMPLES:: 541 542 sage: K.<a> = NumberField(x^3x1) 543 sage: chebyshev_T(5, a) 544 16*a^2 + a  4 545 sage: chebyshev_T(5,MatrixSpace(ZZ, 2)([1, 2, 4, 7])) 546 [40799 44162] 547 [88324 91687] 548 sage: R.<x> = QQ[] 549 sage: parent(chebyshev_T(5, x)) 550 Univariate Polynomial Ring in x over Rational Field 551 """ 552 553 if n in ZZ: # check if n is integer > consider polynomial as algebraic structure 554 return self._eval_(n,x) # Let eval methode decide which is best 555 else: # Consider OrthogonalPolynomial as symbol 556 return super(OrthogonalPolynomial,self).__call__(n,x) 557 547 558 548 559 class Func_chebyshev_T(OrthogonalPolynomial): 549 560 """ … … 579 590 OrthogonalPolynomial.__init__(self, "chebyshev_T", nargs=2, 580 591 conversions=dict(maxima='chebyshev_t', 581 592 mathematica='ChebyshevT')) 582 593 583 594 def _eval_special_values_(self,*args): 584 595 """ 585 596 Values known for special values of x. … … 603 614 if (args[1] == 0): 604 615 return (1+(1)**args[0])*(1)**(args[0]/2)/2 605 616 617 if args[0] < 0 and args[0] in ZZ: 618 return self._eval_(args[0],args[1]) 619 606 620 raise ValueError("value not found") 607 621 608 622 def _evalf_(self, *args,**kwds): … … 617 631 3363.00000000000 618 632 sage: chebyshev_T(5,0.3).n() 619 633 0.998880000000000 620 sage: parent(chebyshev_T(10^2, RIF(2))) # Test for correct data type621 Real Interval Field with 53 bits of precision622 634 """ 623 635 try: 624 636 step_parent = kwds['parent'] … … 648 660 x = args[1] 649 661 return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x}) 650 662 651 def _cheb_recur_(self,n, x, both=False):652 """653 Method for the generalized recursion formula for Chebyshev polynomials,654 suggested by Frederik Johansson.655 returns (T(n,x), T(n1,x)), or (T(n,x), _) if both=False656 657 EXAMPLES::658 659 sage: chebyshev_T._cheb_recur_(5,x)660 (2*(2*(2*x^2  1)*x  x)*(2*x^2  1)  x, False)661 """662 if n == 0:663 return 1, x664 if n == 1:665 return x, 1666 a, b = self._cheb_recur_((n+1)//2, x, both or n % 2)667 if n % 2 == 0:668 return 2*a**2  1, both and 2*a*b  x669 else:670 return 2*a*b  x, both and 2*b**2  1671 672 673 663 def _clenshaw_method_(self,*args): 674 664 """ 675 665 Clenshaw method for :class:`chebyshev_T` (means use recursions in this 676 666 case). This is much faster for numerical evaluation than maxima! 677 See [ASHandbook]_ 227 (p. 782) for details for the recurions 667 See [ASHandbook]_ 227 (p. 782) for details for the recurions. 668 See also [EffCheby] for fast evaluation techniques. 678 669 679 670 EXAMPLES:: 680 671 681 672 sage: chebyshev_T._clenshaw_method_(2,0.1) == chebyshev_T._evalf_(2,0.1) 682 673 True 674 sage: chebyshev_T(51,x) 675 2*(2*(2*(2*(2*(2*x^2  1)^2  1)*(2*(2*x^2  1)*x  x)  x)*(2*(2*(2*x^2  1)*x  x)^2  1)  x)^2  1)*(2*(2*(2*(2*(2*x^2  1)^2  1)*(2*(2*x^2  1)*x  x)  x)*(2*(2*(2*x^2  1)*x  x)^2  1)  x)*(2*(2*(2*(2*x^2  1)*x  x)^2  1)^2  1)  x)  x 676 sage: chebyshev_T._clenshaw_method_(10,x) 677 512*x^10  1280*x^8 + 1120*x^6  400*x^4 + 50*x^2  1 678 683 679 """ 684 n= args[0]680 k = args[0] 685 681 x = args[1] 686 682 687 return self._cheb_recur_(n, x, False)[0] 688 # if k == 0: 689 # return 1 690 # if k == 1: 691 # return x 683 if k == 0: 684 return 1 685 if k == 1: 686 return x 692 687 693 688 # TODO: When evaluation of Symbolic Expressions works better 694 689 # use these explicit formulas instead! … … 698 693 # return cosh(k*acosh(x)) 699 694 #else: # x < 1 700 695 # return (1)**(k%2)*cosh(k*acosh(x)) 696 697 help1 = 1 698 help2 = x 699 if is_Expression(x) and k <= 50: 700 # Recursion gives more compact representations for large k 701 help3 = 0 702 for j in xrange(0,floor(k/2)+1): 703 f = factorial(kj1) / factorial(j) / factorial(k2*j) 704 help3 = help3 + (1)**j * (2*x)**(k2*j) * f 705 help3 = help3 * k / 2 706 return help3 707 else: 708 return self._cheb_recur_(k,x)[0] 709 710 def _cheb_recur_(self,n, x, both=False): 711 """ 712 Generalized recursion formula for Chebyshev polynomials. 713 Implementation suggested by Frederik Johansson. 714 returns (T(n,x), T(n1,x)), or (T(n,x), _) if both=False 715 716 EXAMPLES:: 717 718 sage: chebyshev_T._cheb_recur_(5,x) 719 (2*(2*(2*x^2  1)*x  x)*(2*x^2  1)  x, False) 720 """ 721 if n == 0: 722 return 1, x 723 if n == 1: 724 return x, 1 725 a, b = self._cheb_recur_((n+1)//2, x, both or n % 2) 726 if n % 2 == 0: 727 return 2*a**2  1, both and 2*a*b  x 728 else: 729 return 2*a*b  x, both and 2*b**2  1 701 730 702 731 703 732 def _eval_numpy_(self, *args): … … 711 740 sage: z2 = numpy.array([[1,2],[1,2]]) 712 741 sage: z3 = numpy.array([1,2,3.]) 713 742 sage: chebyshev_T(1,z) 714 array([ 1., 2.])743 array([1, 2]) 715 744 sage: chebyshev_T(1,z2) 716 array([[ 1., 2.],717 [ 1., 2.]])745 array([[1, 2], 746 [1, 2]]) 718 747 sage: chebyshev_T(1,z3) 719 748 array([ 1., 2., 3.]) 720 749 sage: chebyshev_T(z,0.1) … … 739 768 sage: derivative(chebyshev_T(k,x),k) 740 769 Traceback (most recent call last): 741 770 ... 742 ArithmeticError: derivative w.r.t. to the index is not supported yet771 NotImplementedError: derivative w.r.t. to the index is not supported yet 743 772 """ 744 773 diff_param = kwds['diff_param'] 745 774 if diff_param == 0: 746 raise ArithmeticError("derivative w.r.t. to the index is not supported yet")775 raise NotImplementedError("derivative w.r.t. to the index is not supported yet") 747 776 748 777 return args[0]*chebyshev_U(args[0]1,args[1]) 749 778 … … 779 808 OrthogonalPolynomial.__init__(self, "chebyshev_U", nargs=2, 780 809 conversions=dict(maxima='chebyshev_u', 781 810 mathematica='ChebyshevU')) 782 783 def _cheb_recur_(self,n, x, both=False):784 """785 Method for the generalized recursion formula for Chebyshev polynomials,786 suggested by Frederik Johansson.787 returns (T(n,x), T(n1,x)), or (T(n,x), _) if both=False788 789 EXAMPLES::790 791 sage: chebyshev_U._cheb_recur_(3,x)792 (4*(2*x^2  1)*x, False)793 794 """795 796 if n == 0:797 return 1, both and 2*x798 if n == 1:799 return 2*x, both and 4*x**21800 801 a, b = self._cheb_recur_((n1)//2, x, True)802 if n % 2 == 0:803 return (b+a)*(ba), both and 2*b*(x*ba)804 else:805 return 2*a*(bx*a), both and (b+a)*(ba)806 return self._cheb_recur_(n, x, False)[0]807 811 808 812 def _clenshaw_method_(self,*args): 809 813 """ 810 814 Clenshaw method for :class:`chebyshev_U` (means we use the recursion) 811 815 This is much faster for numerical evaluation than maxima. 812 816 See [ASHandbook]_ 227 (p. 782) for details on the recurions. 817 See also [ 813 818 814 819 EXAMPLES:: 815 820 … … 818 823 """ 819 824 k = args[0] 820 825 x = args[1] 826 827 if k == 0: 828 return 1 829 if k == 1: 830 return 2*x 821 831 822 return self._cheb_recur_(k,x)[0] 832 help1 = 1 833 help2 = 2*x 834 if is_Expression(x) and k <= 50: 835 # Recursion gives more compact representations for large k 836 help3 = 0 837 for j in xrange(0,floor(k/2)+1): 838 f = factorial(kj) / factorial(j) / factorial(k2*j) # Change to a binomial? 839 help3 = help3 + (1)**j * (2*x)**(k2*j) * f 840 return help3 841 842 else: 843 return self._cheb_recur_(k,x)[0] 844 845 846 def _cheb_recur_(self,n, x, both=False): 847 """ 848 Generalized recursion formula for Chebyshev polynomials. 849 Implementation suggested by Frederik Johansson. 850 returns (U(n,x), U(n1,x)), or (U(n,x), _) if both=False 851 852 EXAMPLES:: 853 854 sage: chebyshev_U._cheb_recur_(3,x) 855 (4*(2*x^2  1)*x, False) 856 sage: chebyshev_U._cheb_recur_(5,x)[0] 857 2*((2*x + 1)*(2*x  1)*x  4*(2*x^2  1)*x)*(2*x + 1)*(2*x  1) 858 sage: abs(pari('polchebyshev(5, 2, 0.1)')  chebyshev_U(5,0.1)) < 1e10 859 True 860 861 862 """ 863 864 if n == 0: 865 return 1, both and 2*x 866 if n == 1: 867 return 2*x, both and 4*x**21 868 869 a, b = self._cheb_recur_((n1)//2, x, True) 870 if n % 2 == 0: 871 return (b+a)*(ba), both and 2*b*(x*ba) 872 else: 873 return 2*a*(bx*a), both and (b+a)*(ba) 823 874 824 875 def _maxima_init_evaled_(self, *args): 825 876 """ … … 834 885 n = args[0] 835 886 x = args[1] 836 887 return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x}) 837 888 838 889 def _evalf_(self, *args,**kwds): 839 890 """ 840 891 Evaluate :class:`chebyshev_U` numerically with mpmath. … … 853 904 854 905 try: 855 906 precision = step_parent.prec() 856 except ValueError:907 except AttributeError: 857 908 precision = RR.prec() 858 909 859 910 from sage.libs.mpmath.all import call as mpcall … … 883 934 if (args[1] == 0): 884 935 return (1+(1)**args[0])*(1)**(args[0]/2)/2 885 936 937 if args[0] < 0 and args[0] in ZZ: 938 return self._eval_(args[0]2,args[1]) 939 886 940 raise ValueError("value not found") 887 941 888 942 def _eval_numpy_(self, *args): … … 896 950 sage: z2 = numpy.array([[1,2],[1,2]]) 897 951 sage: z3 = numpy.array([1,2,3.]) 898 952 sage: chebyshev_U(1,z) 899 array([ 2., 4.])953 array([2, 4]) 900 954 sage: chebyshev_U(1,z2) 901 array([[ 2., 4.],902 [ 2., 4.]])955 array([[2, 4], 956 [2, 4]]) 903 957 sage: chebyshev_U(1,z3) 904 958 array([ 2., 4., 6.]) 905 959 sage: chebyshev_U(z,0.1) … … 925 979 sage: derivative(chebyshev_U(k,x),k) 926 980 Traceback (most recent call last): 927 981 ... 928 ArithmeticError: derivative w.r.t. to the index is not supported yet982 NotImplementedError: derivative w.r.t. to the index is not supported yet 929 983 """ 930 984 diff_param = kwds['diff_param'] 931 985 if diff_param == 0: 932 raise ArithmeticError("derivative w.r.t. to the index is not supported yet")986 raise NotImplementedError("derivative w.r.t. to the index is not supported yet") 933 987 934 988 return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])args[1]* 935 989 chebyshev_U(args[0],args[1]))/(args[1]**21)