Ticket #9706: trac_9706_chebyshev_corrections.patch

File trac_9706_chebyshev_corrections.patch, 17.3 KB (added by maldun, 7 years ago)

corrected name of file

  • sage/functions/orthogonal_polys.py

    # HG changeset patch
    # User Stefan Reiterer <domors@gmx.net>
    # Date 1386536477 -3600
    # Node ID 76e0bcf8d3bc159779269cb823abf4efd4bbfc33
    # Parent  2419063308f41fb294d04b917d1d1184ed876526
    trac 9706: further improvements, bugfixes and more doctests.
    
    diff --git a/sage/functions/orthogonal_polys.py b/sage/functions/orthogonal_polys.py
    a b  
    315315#                  http://www.gnu.org/licenses/
    316316#*****************************************************************************
    317317
     318import warnings
     319
    318320from sage.misc.sage_eval import sage_eval
    319 from sage.rings.all import ZZ, RR
     321from sage.rings.all import ZZ, RR, CC, RIF, CIF
    320322from sage.calculus.calculus import maxima
    321323
    322324
     
    399401        except KeyError:
    400402            self._maxima_name = None
    401403
    402         BuiltinFunction.__init__(self, name=name, nargs=nargs,
     404        super(OrthogonalPolynomial,self).__init__(name=name, nargs=nargs,
    403405                                 latex_name=latex_name, conversions=conversions)
    404406
    405407    def _maxima_init_evaled_(self, *args):
     
    420422        """
    421423        return None
    422424
    423     def _clenshaw_method_(self, *args):
     425    def _apply_formula_(self, *args):
    424426        """
    425427        The Clenshaw method uses the three term recursion of the polynomial,
    426428        or explicit formulas instead of maxima to evaluate the polynomial
     
    432434
    433435            sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial
    434436            sage: new = OrthogonalPolynomial('testo_P')
    435             sage: new._clenshaw_method_(1,2.0)
     437            sage: new._apply_formula_(1,2.0)
    436438            Traceback (most recent call last):
    437439            ...
    438440            NotImplementedError: no recursive calculation of values implemented
     
    474476            64*x^7 - 112*x^5 + 56*x^3 - 7*x
    475477            sage: chebyshev_T(3/2,x)
    476478            chebyshev_T(3/2, x)
     479            sage: x = PolynomialRing(QQ, 'x').gen()
     480            sage: chebyshev_T(2,x)
     481            2*x^2 - 1
     482            sage: chebyshev_U(2,x)
     483            4*x^2 - 1
    477484            sage: parent(chebyshev_T(4, RIF(5)))
    478485            Real Interval Field with 53 bits of precision
     486            sage: RR2 = RealField(5)
     487            sage: chebyshev_T(100000,RR2(2))
     488            8.9e57180
     489            sage: chebyshev_T(5,Qp(3)(2))
     490            2 + 3^2 + 3^3 + 3^4 + 3^5 + O(3^20)
     491            sage: chebyshev_T(100001/2, 2)
     492            doctest:507: RuntimeWarning: Warning: mpmath returns NoConvergence exception! Use other method instead.
     493            chebyshev_T(100001/2, 2)
     494            sage: chebyshev_U._eval_(1.5, Mod(8,9)) is None
     495            True
    479496        """
    480497        if not is_Expression(args[0]):
    481             if not is_Expression(args[-1]) and is_inexact(args[-1]):
     498        # If x is no expression and is inexact or n is not an integer -> make numerical evaluation
     499            if (not is_Expression(args[-1])) and (is_inexact(args[-1]) or not args[0] in ZZ):
    482500                try:
    483501                    import sage.libs.mpmath.all as mpmath
    484                     step_parent = parent(args[-1])
    485                     return step_parent(self._evalf_(*args))
     502                    return self._evalf_(*args)
    486503                except AttributeError:
    487504                    pass
    488505                except mpmath.NoConvergence:
    489                     print "Warning: mpmath returns NoConvergence!"
    490                     print "Switching to clenshaw_method, but it may not be stable!"
     506                    warnings.warn("Warning: mpmath returns NoConvergence exception! Use other method instead.",
     507                                  RuntimeWarning)
    491508                except ValueError:
    492509                    pass
    493510
    494             # A faster check would be nice...
    495             if args[0] != floor(args[0]):
    496                 if not is_Expression(args[-1]):
    497                     try:
    498                         return self._evalf_(*args)
    499                     except AttributeError:
    500                         pass
    501                 else:
     511            # n is not an integer and x is an expression -> return symbolic expression.
     512            if not args[0] in ZZ:
     513                if is_Expression(args[-1]):
    502514                    return None
     515
     516        # Check for known identities
    503517        try:
    504518            return self._eval_special_values_(*args)
    505519        except ValueError:
     
    511525        if args[0] < 0 and args[0] in ZZ:
    512526                return None
    513527
    514         if not is_Expression(args[0]):
     528        if args[0] in ZZ:
    515529            try:
    516                 return self._clenshaw_method_(*args)
     530                return self._apply_formula_(*args)
    517531            except NotImplementedError:
    518532                pass
    519533
    520534        if self._maxima_name is None:
    521535            return None
    522536
    523         _init()
    524         try: # use maxima as last resort
    525             s = maxima(self._maxima_init_evaled_(*args))
    526             return self._maxima_init_evaled_(*args)
    527         except TypeError:
     537        if args[0] in ZZ: # use maxima as last resort
     538            return self._old_maxima_(*args)
     539        else:
    528540            return None
    529541       
    530         if self._maxima_name in repr(s):
    531             return None
    532 
    533         return s.sage()
    534 
    535     def __call__(self,n,x):
     542    def __call__(self,*args,**kwds):
    536543        """
    537544        This overides the call method from SageObject to avoid problems with coercions,
    538         since the _eval_ method is able to handle for more data types than symbolic functions
     545        since the _eval_ method is able to handle more data types than symbolic functions
    539546        would normally allow.
    540547        Thus we have the distinction between algebraic objects (if n is an integer),
    541548        and else as symbolic function.
     
    552559            sage: parent(chebyshev_T(5, x))
    553560            Univariate Polynomial Ring in x over Rational Field
    554561        """
     562        if 'hold' not in kwds:
     563            kwds['hold'] = False
     564        if 'coerce' not in kwds:
     565            kwds['coerce']=True
    555566       
    556         if n in ZZ: # check if n is integer -> consider polynomial as algebraic structure
    557             return self._eval_(n,x) # Let eval methode decide which is best
     567        if args[0] in ZZ and kwds['hold'] is False: #check if n is in ZZ->consider polynomial as algebraic structure
     568            return self._eval_(*args) # Let eval methode decide which is best
    558569        else: # Consider OrthogonalPolynomial as symbol
    559             return super(OrthogonalPolynomial,self).__call__(n,x)
     570            return super(OrthogonalPolynomial,self).__call__(*args,**kwds)
    560571
     572    def _old_maxima_(self,*args):
     573        """
     574        Method which holds the old maxima wrappers as last alternative.
     575        It returns None per default, and it only needs to be implemented,
     576        if it is necessary.
     577
     578        EXAMPLES::
     579       
     580            sage: chebyshev_T._old_maxima_(-7,x) is None
     581            True
     582        """
     583        None
    561584
    562585class Func_chebyshev_T(OrthogonalPolynomial):
    563586    """
     
    590613            sage: chebyshev_T2(1,x)
    591614            x
    592615        """
    593         OrthogonalPolynomial.__init__(self, "chebyshev_T", nargs=2,
     616        super(Func_chebyshev_T,self).__init__("chebyshev_T", nargs=2,
    594617                                      conversions=dict(maxima='chebyshev_t',
    595618                                                       mathematica='ChebyshevT'))
    596619   
     
    607630            1
    608631            sage: chebyshev_T(n,-1)
    609632            (-1)^n
     633            sage: chebyshev_T(-7, x) - chebyshev_T(7,x)
     634            0
     635            sage: chebyshev_T._eval_special_values_(3/2,x)
     636            Traceback (most recent call last):
     637            ...
     638            ValueError: No special values for non integral indices!
     639            sage: chebyshev_T._eval_special_values_(n, 0.1)
     640            Traceback (most recent call last):
     641            ...
     642            ValueError: Value not found!
     643            sage: chebyshev_T._eval_special_values_(26, Mod(9,9))
     644            Traceback (most recent call last):
     645            ...
     646            ValueError: Value not found!
    610647        """
     648        if (not is_Expression(args[0])) and (not args[0] in ZZ):
     649            raise ValueError("No special values for non integral indices!")
     650       
    611651        if args[-1] == 1:
    612             return 1
     652            return args[-1]
    613653       
    614654        if args[-1] == -1:
    615             return (-1)**args[0]
     655            return args[-1]**args[0]
    616656
    617         if (args[-1] == 0):
     657        if (args[-1] == 0 and args[-1] in CC):
    618658            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
    619659
    620660        if args[0] < 0 and args[0] in ZZ:
    621661            return self._eval_(-args[0],args[-1])
    622662
    623         raise ValueError("value not found")
     663        raise ValueError("Value not found!")
    624664   
    625665    def _evalf_(self, *args,**kwds):
    626666        """
    627667        Evaluates :class:`chebyshev_T` numerically with mpmath.
     668        If the index is an integer we use the recursive formula since
     669        it is faster.
    628670
    629671        EXAMPLES::
    630672
     
    636678            0.998880000000000
    637679            sage: chebyshev_T(1/2, 0)
    638680            0.707106781186548
     681            sage: chebyshev_T._evalf_(1.5, Mod(8,9))
     682            Traceback (most recent call last):
     683            ...
     684            ValueError: No compatible type!
     685
    639686        """
     687        if args[0] in ZZ:
     688            return self._cheb_recur_(*args)[0]
     689       
    640690        try:
    641             step_parent = kwds['parent']
     691            real_parent = kwds['parent']
    642692        except KeyError:
    643             step_parent = parent(args[-1])
     693            real_parent = parent(args[-1])
    644694
    645         try:
    646             precision = step_parent.prec()
    647         except AttributeError:
    648             precision = RR.prec()
     695        x_set = False
     696        if hasattr(real_parent,"precision"): # Check if we have a data type with precision
     697            x = args[-1]
     698            step_parent = real_parent
     699            x_set = True
     700        else:
     701            if args[-1] in RR:
     702                x = RR(args[-1])
     703                step_parent = RR
     704                x_set = True
     705            elif args[-1] in CC:
     706                x = CC(args[-1])
     707                step_parent = CC
     708                x_set = True
    649709
     710        if not x_set:
     711            raise ValueError("No compatible type!")
     712               
    650713        from sage.libs.mpmath.all import call as mpcall
    651714        from sage.libs.mpmath.all import chebyt as mpchebyt
    652715
    653         return mpcall(mpchebyt,args[0],args[-1],prec = precision)
     716        return mpcall(mpchebyt,args[0],x,parent=step_parent)
    654717
    655718    def _maxima_init_evaled_(self, *args):
    656719        """
     
    659722        EXAMPLES::
    660723
    661724            sage: chebyshev_T._maxima_init_evaled_(1,x)
    662             x
     725            'x'
     726            sage: var('n')
     727            n
     728            sage: maxima(chebyshev_T(n,x))
     729            chebyshev_t(n,x)
     730
    663731        """
    664732        n = args[0]
    665733        x = args[1]
    666         return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
     734        return maxima.eval('chebyshev_t({0},{1})'.format(n,x))
    667735
    668     def _clenshaw_method_(self,*args):
     736       
     737    def _apply_formula_(self,*args):
    669738        """
    670739        Clenshaw method for :class:`chebyshev_T` (means use recursions in this
    671740        case). This is much faster for numerical evaluation than maxima!
     
    674743
    675744        EXAMPLES::
    676745
    677             sage: chebyshev_T._clenshaw_method_(2,0.1) == chebyshev_T._evalf_(2,0.1)
     746            sage: chebyshev_T._apply_formula_(2,0.1) == chebyshev_T._evalf_(2,0.1)
    678747            True
    679748            sage: chebyshev_T(51,x)
    680749            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
    681             sage: chebyshev_T._clenshaw_method_(10,x)
     750            sage: chebyshev_T._apply_formula_(10,x)
    682751            512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1
    683752
    684753        """
     
    701770       
    702771        help1 = 1
    703772        help2 = x
    704         if is_Expression(x) and k <= 50:
     773        if is_Expression(x) and k <= 25:
    705774            # Recursion gives more compact representations for large k
    706775            help3 = 0
    707776            for j in xrange(0,floor(k/2)+1):
     
    814883                                      conversions=dict(maxima='chebyshev_u',
    815884                                                       mathematica='ChebyshevU'))
    816885       
    817     def _clenshaw_method_(self,*args):
     886    def _apply_formula_(self,*args):
    818887        """
    819888        Clenshaw method for :class:`chebyshev_U` (means we use the recursion)
    820889        This is much faster for numerical evaluation than maxima.
     
    823892
    824893        EXAMPLES::
    825894
    826             sage: chebyshev_U._clenshaw_method_(2,0.1) == chebyshev_U._evalf_(2,0.1)
     895            sage: chebyshev_U._apply_formula_(2,0.1) == chebyshev_U._evalf_(2,0.1)
    827896            True
    828897        """
    829898        k = args[0]
     
    836905
    837906        help1 = 1
    838907        help2 = 2*x
    839         if is_Expression(x) and k <= 50:
     908        if is_Expression(x) and k <= 25:
    840909            # Recursion gives more compact representations for large k
    841910            help3 = 0
    842911            for j in xrange(0,floor(k/2)+1):
     
    862931            -2*((2*x + 1)*(2*x - 1)*x - 4*(2*x^2 - 1)*x)*(2*x + 1)*(2*x - 1)
    863932            sage: abs(pari('polchebyshev(5, 2, 0.1)') - chebyshev_U(5,0.1)) < 1e-10
    864933            True
    865 
    866 
    867934        """
    868935       
    869936        if n == 0:
     
    883950
    884951        EXAMPLES::
    885952
    886             sage: x = PolynomialRing(QQ, 'x').gen()
    887             sage: chebyshev_T(2,x)
    888             2*x^2 - 1
     953            sage: maxima(chebyshev_U(5,x))
     954            32*x^5-32*x^3+6*x
     955            sage: var('n')
     956            n
     957            sage: maxima(chebyshev_U(n,x))
     958            chebyshev_u(n,x)
     959            sage: maxima(chebyshev_U(2,x))
     960            4*x^2-1
    889961        """
    890962        n = args[0]
    891963        x = args[1]
    892         return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
     964        return maxima.eval('chebyshev_u({0},{1})'.format(n,x))
    893965
    894966    def _evalf_(self, *args,**kwds):
    895967        """
    896968        Evaluate :class:`chebyshev_U` numerically with mpmath.
     969        If index is an integer use recursive formula since it is faster,
     970        for chebyshev polynomials.
    897971
    898972        EXAMPLES::
    899973
     
    901975            98280.0000000000 - 11310.0000000000*I
    902976            sage: chebyshev_U(10,3).n(75)
    903977            4.661117900000000000000e7
     978            sage: chebyshev_U._evalf_(1.5, Mod(8,9))
     979            Traceback (most recent call last):
     980            ...
     981            ValueError: No compatible type!
    904982        """
     983        if args[0] in ZZ:
     984            return self._cheb_recur_(*args)[0]
    905985        try:
    906             step_parent = kwds['parent']
     986            real_parent = kwds['parent']
    907987        except KeyError:
    908             step_parent = parent(args[-1])
     988            real_parent = parent(args[-1])
    909989
    910         try:
    911             precision = step_parent.prec()
    912         except AttributeError:
    913             precision = RR.prec()
     990        x_set = False
     991        if hasattr(real_parent,"precision"): # Check if we have a data type with precision
     992            x = args[-1]
     993            step_parent = real_parent
     994            x_set = True
     995        else:
     996            if args[-1] in RR:
     997                x = RR(args[-1])
     998                step_parent = RR
     999                x_set = True
     1000            elif args[-1] in CC:
     1001                x = CC(args[-1])
     1002                step_parent = CC
     1003                x_set = True
     1004
     1005        if not x_set:
     1006            raise ValueError("No compatible type!")
    9141007
    9151008        from sage.libs.mpmath.all import call as mpcall
    9161009        from sage.libs.mpmath.all import chebyu as mpchebyu
    9171010
    918         return mpcall(mpchebyu,args[0],args[-1],prec = precision)
     1011        return mpcall(mpchebyu,args[0],args[-1],parent = step_parent)
    9191012
    9201013    def _eval_special_values_(self,*args):
    9211014        """
     
    9291022            n + 1
    9301023            sage: chebyshev_U(n,-1)
    9311024            (-1)^n*(n + 1)
     1025            sage: chebyshev_U._eval_special_values_(26, Mod(0,9))
     1026            Traceback (most recent call last):
     1027            ...
     1028            ValueError: Value not found!
     1029            sage: parent(chebyshev_U(3, Mod(8,9)))
     1030            Ring of integers modulo 9
     1031            sage: parent(chebyshev_U(3, Mod(1,9)))
     1032            Ring of integers modulo 9
     1033            sage: chebyshev_U(n, 0)
     1034            1/2*(-1)^(1/2*n)*((-1)^n + 1)
     1035            sage: chebyshev_U(-3,x) + chebyshev_U(1,x)
     1036            0
     1037            sage: chebyshev_U._eval_special_values_(1.5, Mod(8,9))
     1038            Traceback (most recent call last):
     1039            ...
     1040            ValueError: No special values for non integral indices!
    9321041        """
     1042        if (not is_Expression(args[0])) and (not args[0] in ZZ):
     1043            raise ValueError("No special values for non integral indices!")
     1044 
    9331045        if args[-1] == 1:
    934             return (args[0]+1)
     1046            return args[-1]*(args[0]+1)
    9351047       
    9361048        if args[-1] == -1:
    937             return (-1)**args[0]*(args[0]+1)
     1049            return args[-1]**args[0]*(args[0]+1)
    9381050
    939         if (args[-1] == 0):
     1051        if (args[-1] == 0 and args[-1] in CC):
    9401052            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
    9411053
    9421054        if args[0] < 0 and args[0] in ZZ:
    9431055            return -self._eval_(-args[0]-2,args[-1])
    9441056
    945         raise ValueError("value not found")
     1057        raise ValueError("Value not found!")
    9461058
    9471059    def _eval_numpy_(self, *args):
    9481060        """