Ticket #9706: trac_9706_chebyshev_patch_new.patch

File trac_9706_chebyshev_patch_new.patch, 14.5 KB (added by maldun, 8 years ago)

Patch in a whole with corrections etc.

  • 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  
    287287
    288288-  :wikipedia:`Associated_Legendre_polynomials`
    289289
     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): 79-99.
    290294
    291295AUTHORS:
    292296
     
    448452        """
    449453        raise ValueError("no special values known")
    450454
    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 functions
    455         would normally allow.
    456 
    457         EXAMPLES::
    458 
    459            sage: K.<a> = NumberField(x^3-x-1)
    460            sage: chebyshev_T(5, a)
    461            16*a^2 + a - 4
    462            
    463          
    464         """
    465         try:
    466             return super(OrthogonalPolynomial,self).__call__(n,x)
    467         except TypeError:
    468              return self._eval_(n,x)
    469 
    470    
    471455    def _eval_(self, *args):
    472456        """
    473457        The symbolic evaluation is either done with maxima, or with direct
     
    487471            sage: chebyshev_T(n,-1)
    488472            (-1)^n
    489473            sage: chebyshev_T(-7,x)
    490             chebyshev_T(-7, x)
     474            64*x^7 - 112*x^5 + 56*x^3 - 7*x
    491475            sage: chebyshev_T(3/2,x)
    492476            chebyshev_T(3/2, x)
    493477        """
     
    495479            if not is_Expression(args[-1]) and is_inexact(args[-1]):
    496480                try:
    497481                    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)
    500483                except AttributeError:
    501484                    pass
    502485                except mpmath.NoConvergence:
     
    514497                        pass
    515498                else:
    516499                    return None
    517        
    518             if args[0] < 0:
    519                 return None
    520            
    521500        try:
    522501            return self._eval_special_values_(*args)
    523502        except ValueError:
    524503            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
    526511        if not is_Expression(args[0]):
    527512            try:
    528513                return self._clenshaw_method_(*args)
     
    544529
    545530        return s.sage()
    546531
     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^3-x-1)
     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
    547558
    548559class Func_chebyshev_T(OrthogonalPolynomial):
    549560    """
     
    579590        OrthogonalPolynomial.__init__(self, "chebyshev_T", nargs=2,
    580591                                      conversions=dict(maxima='chebyshev_t',
    581592                                                       mathematica='ChebyshevT'))
    582        
     593   
    583594    def _eval_special_values_(self,*args):
    584595        """
    585596        Values known for special values of x.
     
    603614        if (args[-1] == 0):
    604615            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
    605616
     617        if args[0] < 0 and args[0] in ZZ:
     618            return self._eval_(-args[0],args[-1])
     619
    606620        raise ValueError("value not found")
    607621   
    608622    def _evalf_(self, *args,**kwds):
     
    617631            -3363.00000000000
    618632            sage: chebyshev_T(5,0.3).n()
    619633            0.998880000000000
    620             sage: parent(chebyshev_T(10^2, RIF(2))) # Test for correct data type
    621             Real Interval Field with 53 bits of precision
    622634        """
    623635        try:
    624636            step_parent = kwds['parent']
     
    648660        x = args[1]
    649661        return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
    650662
    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(n-1,x)), or (T(n,x), _) if both=False
    656 
    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, x
    664         if n == 1:
    665             return x, 1
    666         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 - x
    669         else:
    670             return 2*a*b - x, both and 2*b**2 - 1
    671 
    672    
    673663    def _clenshaw_method_(self,*args):
    674664        """
    675665        Clenshaw method for :class:`chebyshev_T` (means use recursions in this
    676666        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.
    678669
    679670        EXAMPLES::
    680671
    681672            sage: chebyshev_T._clenshaw_method_(2,0.1) == chebyshev_T._evalf_(2,0.1)
    682673            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
    683679        """
    684         n = args[0]
     680        k = args[0]
    685681        x = args[1]
    686682
    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
    692687
    693688        # TODO: When evaluation of Symbolic Expressions works better
    694689        # use these explicit formulas instead!
     
    698693        #    return cosh(k*acosh(x))
    699694        #else: # x < -1
    700695        #    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(k-j-1) / factorial(j) / factorial(k-2*j)
     704                help3 = help3 + (-1)**j * (2*x)**(k-2*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(n-1,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
    701730
    702731
    703732    def _eval_numpy_(self, *args):
     
    711740            sage: z2 = numpy.array([[1,2],[1,2]])
    712741            sage: z3 = numpy.array([1,2,3.])
    713742            sage: chebyshev_T(1,z)
    714             array([ 1.,  2.])
     743            array([1, 2])
    715744            sage: chebyshev_T(1,z2)
    716             array([[ 1.,  2.],
    717                    [ 1.,  2.]])
     745            array([[1, 2],
     746                   [1, 2]])
    718747            sage: chebyshev_T(1,z3)
    719748            array([ 1.,  2.,  3.])
    720749            sage: chebyshev_T(z,0.1)
     
    739768            sage: derivative(chebyshev_T(k,x),k)
    740769            Traceback (most recent call last):
    741770            ...
    742             ArithmeticError: derivative w.r.t. to the index is not supported yet
     771            NotImplementedError: derivative w.r.t. to the index is not supported yet
    743772        """
    744773        diff_param = kwds['diff_param']
    745774        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")
    747776
    748777        return args[0]*chebyshev_U(args[0]-1,args[1])
    749778   
     
    779808        OrthogonalPolynomial.__init__(self, "chebyshev_U", nargs=2,
    780809                                      conversions=dict(maxima='chebyshev_u',
    781810                                                       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(n-1,x)), or (T(n,x), _) if both=False
    788 
    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*x
    798         if n == 1:
    799             return 2*x, both and 4*x**2-1
    800 
    801         a, b = self._cheb_recur_((n-1)//2, x, True)
    802         if n % 2 == 0:
    803             return (b+a)*(b-a), both and 2*b*(x*b-a)
    804         else:
    805             return 2*a*(b-x*a), both and (b+a)*(b-a)
    806         return self._cheb_recur_(n, x, False)[0]
    807811       
    808812    def _clenshaw_method_(self,*args):
    809813        """
    810814        Clenshaw method for :class:`chebyshev_U` (means we use the recursion)
    811815        This is much faster for numerical evaluation than maxima.
    812816        See [ASHandbook]_ 227 (p. 782) for details on the recurions.
     817        See also [
    813818
    814819        EXAMPLES::
    815820
     
    818823        """
    819824        k = args[0]
    820825        x = args[1]
     826           
     827        if k == 0:
     828            return 1
     829        if k == 1:
     830            return 2*x
    821831
    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(k-j) / factorial(j) / factorial(k-2*j) # Change to a binomial?
     839                help3 = help3 + (-1)**j * (2*x)**(k-2*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(n-1,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)) < 1e-10
     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**2-1
     868           
     869        a, b = self._cheb_recur_((n-1)//2, x, True)
     870        if n % 2 == 0:
     871            return (b+a)*(b-a), both and 2*b*(x*b-a)
     872        else:
     873            return 2*a*(b-x*a), both and (b+a)*(b-a)
    823874
    824875    def _maxima_init_evaled_(self, *args):
    825876        """
     
    834885        n = args[0]
    835886        x = args[1]
    836887        return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
    837        
     888
    838889    def _evalf_(self, *args,**kwds):
    839890        """
    840891        Evaluate :class:`chebyshev_U` numerically with mpmath.
     
    853904
    854905        try:
    855906            precision = step_parent.prec()
    856         except ValueError:
     907        except AttributeError:
    857908            precision = RR.prec()
    858909
    859910        from sage.libs.mpmath.all import call as mpcall
     
    883934        if (args[-1] == 0):
    884935            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
    885936
     937        if args[0] < 0 and args[0] in ZZ:
     938            return -self._eval_(-args[0]-2,args[-1])
     939
    886940        raise ValueError("value not found")
    887941
    888942    def _eval_numpy_(self, *args):
     
    896950            sage: z2 = numpy.array([[1,2],[1,2]])
    897951            sage: z3 = numpy.array([1,2,3.])
    898952            sage: chebyshev_U(1,z)
    899             array([ 2.,  4.])
     953            array([2, 4])
    900954            sage: chebyshev_U(1,z2)
    901             array([[ 2.,  4.],
    902                    [ 2.,  4.]])
     955            array([[2, 4],
     956                   [2, 4]])
    903957            sage: chebyshev_U(1,z3)
    904958            array([ 2.,  4.,  6.])
    905959            sage: chebyshev_U(z,0.1)
     
    925979            sage: derivative(chebyshev_U(k,x),k)
    926980            Traceback (most recent call last):
    927981            ...
    928             ArithmeticError: derivative w.r.t. to the index is not supported yet
     982            NotImplementedError: derivative w.r.t. to the index is not supported yet
    929983        """
    930984        diff_param = kwds['diff_param']
    931985        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")
    933987
    934988        return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
    935989                chebyshev_U(args[0],args[1]))/(args[1]**2-1)