Ticket #9706: trac_9706_new_clenshaw.patch

File trac_9706_new_clenshaw.patch, 7.9 KB (added by maldun, 8 years ago)

new fixes and patches

  • sage/functions/orthogonal_polys.py

    # HG changeset patch
    # User Stefan Reiterer <domors@gmx.net>
    # Date 1386199660 -3600
    # Node ID 12e04e4ed442d2a6b6d9af7919f87ee49191776d
    # Parent  7e924548b0a7891ccec1672e4b7b75d7ccd54ba4
    Trac 9706: better clenshaw, bugfixes
    
    diff --git a/sage/functions/orthogonal_polys.py b/sage/functions/orthogonal_polys.py
    a b  
    448448        """
    449449        raise ValueError("no special values known")
    450450
     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   
    451471    def _eval_(self, *args):
    452472        """
    453473        The symbolic evaluation is either done with maxima, or with direct
     
    475495            if not is_Expression(args[-1]) and is_inexact(args[-1]):
    476496                try:
    477497                    import sage.libs.mpmath.all as mpmath
    478                     return self._evalf_(*args)
     498                    father = parent(args[-1]) # convert back to original type
     499                    return father(self._evalf_(*args))
    479500                except AttributeError:
    480501                    pass
    481502                except mpmath.NoConvergence:
     
    558579        OrthogonalPolynomial.__init__(self, "chebyshev_T", nargs=2,
    559580                                      conversions=dict(maxima='chebyshev_t',
    560581                                                       mathematica='ChebyshevT'))
    561    
     582       
    562583    def _eval_special_values_(self,*args):
    563584        """
    564585        Values known for special values of x.
     
    596617            -3363.00000000000
    597618            sage: chebyshev_T(5,0.3).n()
    598619            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
    599622        """
    600623        try:
    601624            step_parent = kwds['parent']
     
    625648        x = args[1]
    626649        return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
    627650
     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   
    628673    def _clenshaw_method_(self,*args):
    629674        """
    630675        Clenshaw method for :class:`chebyshev_T` (means use recursions in this
     
    636681            sage: chebyshev_T._clenshaw_method_(2,0.1) == chebyshev_T._evalf_(2,0.1)
    637682            True
    638683        """
    639         k = args[0]
     684        n = args[0]
    640685        x = args[1]
    641686
    642         if k == 0:
    643             return 1
    644         if k == 1:
    645             return x
     687        return self._cheb_recur_(n, x, False)[0]
     688        # if k == 0:
     689        #     return 1
     690        # if k == 1:
     691        #     return x
    646692
    647693        # TODO: When evaluation of Symbolic Expressions works better
    648694        # use these explicit formulas instead!
     
    652698        #    return cosh(k*acosh(x))
    653699        #else: # x < -1
    654700        #    return (-1)**(k%2)*cosh(k*acosh(-x))
    655        
    656         help1 = 1
    657         help2 = x
    658         if is_Expression(x):
    659             #raise NotImplementedError
    660             help3 = 0
    661             for j in xrange(0,floor(k/2)+1):
    662                 f = factorial(k-j-1) / factorial(j) / factorial(k-2*j)
    663                 help3 = help3 + (-1)**j * (2*x)**(k-2*j) * f
    664             help3 = help3 * k / 2
    665         else:
    666             for j in xrange(0,k-1):
    667                 help3 = 2*x*help2 - help1
    668                 help1 = help2
    669                 help2 = help3
    670        
    671         return help3
     701
    672702
    673703    def _eval_numpy_(self, *args):
    674704        """
     
    709739            sage: derivative(chebyshev_T(k,x),k)
    710740            Traceback (most recent call last):
    711741            ...
    712             NotImplementedError: derivative w.r.t. to the index is not supported yet
     742            ArithmeticError: derivative w.r.t. to the index is not supported yet
    713743        """
    714744        diff_param = kwds['diff_param']
    715745        if diff_param == 0:
    716             raise NotImplementedError("derivative w.r.t. to the index is not supported yet")
     746            raise ArithmeticError("derivative w.r.t. to the index is not supported yet")
    717747
    718748        return args[0]*chebyshev_U(args[0]-1,args[1])
    719749   
     
    749779        OrthogonalPolynomial.__init__(self, "chebyshev_U", nargs=2,
    750780                                      conversions=dict(maxima='chebyshev_u',
    751781                                                       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]
    752807       
    753808    def _clenshaw_method_(self,*args):
    754809        """
     
    763818        """
    764819        k = args[0]
    765820        x = args[1]
    766            
    767         if k == 0:
    768             return 1
    769         if k == 1:
    770             return 2*x
    771821
    772         help1 = 1
    773         help2 = 2*x
    774         if is_Expression(x):
    775             #raise NotImplementedError("Maxima is faster here!")
    776             help3 = 0
    777             for j in xrange(0,floor(k/2)+1):
    778                 f = factorial(k-j) / factorial(j) / factorial(k-2*j) # Change to a binomial?
    779                 help3 = help3 + (-1)**j * (2*x)**(k-2*j) * f
    780            
    781         else:
    782             for j in xrange(0,k-1):
    783                 help3 = 2*x*help2 - help1
    784                 help1 = help2
    785                 help2 = help3
    786                
    787         return help3   
     822        return self._cheb_recur_(k,x)[0]
    788823
    789824    def _maxima_init_evaled_(self, *args):
    790825        """
     
    818853
    819854        try:
    820855            precision = step_parent.prec()
    821         except AttributeError:
     856        except ValueError:
    822857            precision = RR.prec()
    823858
    824859        from sage.libs.mpmath.all import call as mpcall
     
    890925            sage: derivative(chebyshev_U(k,x),k)
    891926            Traceback (most recent call last):
    892927            ...
    893             NotImplementedError: derivative w.r.t. to the index is not supported yet
     928            ArithmeticError: derivative w.r.t. to the index is not supported yet
    894929        """
    895930        diff_param = kwds['diff_param']
    896931        if diff_param == 0:
    897             raise NotImplementedError("derivative w.r.t. to the index is not supported yet")
     932            raise ArithmeticError("derivative w.r.t. to the index is not supported yet")
    898933
    899934        return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
    900935                chebyshev_U(args[0],args[1]))/(args[1]**2-1)