# Ticket #9706: trac_9707_chebyshev_corrections.patch

File trac_9707_chebyshev_corrections.patch, 17.3 KB (added by maldun, 8 years ago)

Bunch of corrections and new doctests

• ## 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 #                  http://www.gnu.org/licenses/ #***************************************************************************** import warnings from sage.misc.sage_eval import sage_eval from sage.rings.all import ZZ, RR from sage.rings.all import ZZ, RR, CC, RIF, CIF from sage.calculus.calculus import maxima except KeyError: self._maxima_name = None BuiltinFunction.__init__(self, name=name, nargs=nargs, super(OrthogonalPolynomial,self).__init__(name=name, nargs=nargs, latex_name=latex_name, conversions=conversions) def _maxima_init_evaled_(self, *args): """ return None def _clenshaw_method_(self, *args): def _apply_formula_(self, *args): """ The Clenshaw method uses the three term recursion of the polynomial, or explicit formulas instead of maxima to evaluate the polynomial sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial sage: new = OrthogonalPolynomial('testo_P') sage: new._clenshaw_method_(1,2.0) sage: new._apply_formula_(1,2.0) Traceback (most recent call last): ... NotImplementedError: no recursive calculation of values implemented 64*x^7 - 112*x^5 + 56*x^3 - 7*x sage: chebyshev_T(3/2,x) chebyshev_T(3/2, x) sage: x = PolynomialRing(QQ, 'x').gen() sage: chebyshev_T(2,x) 2*x^2 - 1 sage: chebyshev_U(2,x) 4*x^2 - 1 sage: parent(chebyshev_T(4, RIF(5))) Real Interval Field with 53 bits of precision sage: RR2 = RealField(5) sage: chebyshev_T(100000,RR2(2)) 8.9e57180 sage: chebyshev_T(5,Qp(3)(2)) 2 + 3^2 + 3^3 + 3^4 + 3^5 + O(3^20) sage: chebyshev_T(100001/2, 2) doctest:507: RuntimeWarning: Warning: mpmath returns NoConvergence exception! Use other method instead. chebyshev_T(100001/2, 2) sage: chebyshev_U._eval_(1.5, Mod(8,9)) is None True """ if not is_Expression(args[0]): if not is_Expression(args[-1]) and is_inexact(args[-1]): # If x is no expression and is inexact or n is not an integer -> make numerical evaluation if (not is_Expression(args[-1])) and (is_inexact(args[-1]) or not args[0] in ZZ): try: import sage.libs.mpmath.all as mpmath step_parent = parent(args[-1]) return step_parent(self._evalf_(*args)) return self._evalf_(*args) except AttributeError: pass except mpmath.NoConvergence: print "Warning: mpmath returns NoConvergence!" print "Switching to clenshaw_method, but it may not be stable!" warnings.warn("Warning: mpmath returns NoConvergence exception! Use other method instead.", RuntimeWarning) except ValueError: pass # A faster check would be nice... if args[0] != floor(args[0]): if not is_Expression(args[-1]): try: return self._evalf_(*args) except AttributeError: pass else: # n is not an integer and x is an expression -> return symbolic expression. if not args[0] in ZZ: if is_Expression(args[-1]): return None # Check for known identities try: return self._eval_special_values_(*args) except ValueError: if args[0] < 0 and args[0] in ZZ: return None if not is_Expression(args[0]): if args[0] in ZZ: try: return self._clenshaw_method_(*args) return self._apply_formula_(*args) except NotImplementedError: pass if self._maxima_name is None: return None _init() try: # use maxima as last resort s = maxima(self._maxima_init_evaled_(*args)) return self._maxima_init_evaled_(*args) except TypeError: if args[0] in ZZ: # use maxima as last resort return self._old_maxima_(*args) else: return None if self._maxima_name in repr(s): return None return s.sage() def __call__(self,n,x): def __call__(self,*args,**kwds): """ This overides the call method from SageObject to avoid problems with coercions, since the _eval_ method is able to handle for more data types than symbolic functions since the _eval_ method is able to handle more data types than symbolic functions would normally allow. Thus we have the distinction between algebraic objects (if n is an integer), and else as symbolic function. sage: parent(chebyshev_T(5, x)) Univariate Polynomial Ring in x over Rational Field """ if 'hold' not in kwds: kwds['hold'] = False if 'coerce' not in kwds: kwds['coerce']=True if n in ZZ: # check if n is integer -> consider polynomial as algebraic structure return self._eval_(n,x) # Let eval methode decide which is best if args[0] in ZZ and kwds['hold'] is False: #check if n is in ZZ->consider polynomial as algebraic structure return self._eval_(*args) # Let eval methode decide which is best else: # Consider OrthogonalPolynomial as symbol return super(OrthogonalPolynomial,self).__call__(n,x) return super(OrthogonalPolynomial,self).__call__(*args,**kwds) def _old_maxima_(self,*args): """ Method which holds the old maxima wrappers as last alternative. It returns None per default, and it only needs to be implemented, if it is necessary. EXAMPLES:: sage: chebyshev_T._old_maxima_(-7,x) is None True """ None class Func_chebyshev_T(OrthogonalPolynomial): """ sage: chebyshev_T2(1,x) x """ OrthogonalPolynomial.__init__(self, "chebyshev_T", nargs=2, super(Func_chebyshev_T,self).__init__("chebyshev_T", nargs=2, conversions=dict(maxima='chebyshev_t', mathematica='ChebyshevT')) 1 sage: chebyshev_T(n,-1) (-1)^n sage: chebyshev_T(-7, x) - chebyshev_T(7,x) 0 sage: chebyshev_T._eval_special_values_(3/2,x) Traceback (most recent call last): ... ValueError: No special values for non integral indices! sage: chebyshev_T._eval_special_values_(n, 0.1) Traceback (most recent call last): ... ValueError: Value not found! sage: chebyshev_T._eval_special_values_(26, Mod(9,9)) Traceback (most recent call last): ... ValueError: Value not found! """ if (not is_Expression(args[0])) and (not args[0] in ZZ): raise ValueError("No special values for non integral indices!") if args[-1] == 1: return 1 return args[-1] if args[-1] == -1: return (-1)**args[0] return args[-1]**args[0] if (args[-1] == 0): if (args[-1] == 0 and args[-1] in CC): return (1+(-1)**args[0])*(-1)**(args[0]/2)/2 if args[0] < 0 and args[0] in ZZ: return self._eval_(-args[0],args[-1]) raise ValueError("value not found") raise ValueError("Value not found!") def _evalf_(self, *args,**kwds): """ Evaluates :class:`chebyshev_T` numerically with mpmath. If the index is an integer we use the recursive formula since it is faster. EXAMPLES:: 0.998880000000000 sage: chebyshev_T(1/2, 0) 0.707106781186548 sage: chebyshev_T._evalf_(1.5, Mod(8,9)) Traceback (most recent call last): ... ValueError: No compatible type! """ if args[0] in ZZ: return self._cheb_recur_(*args)[0] try: step_parent = kwds['parent'] real_parent = kwds['parent'] except KeyError: step_parent = parent(args[-1]) real_parent = parent(args[-1]) try: precision = step_parent.prec() except AttributeError: precision = RR.prec() x_set = False if hasattr(real_parent,"precision"): # Check if we have a data type with precision x = args[-1] step_parent = real_parent x_set = True else: if args[-1] in RR: x = RR(args[-1]) step_parent = RR x_set = True elif args[-1] in CC: x = CC(args[-1]) step_parent = CC x_set = True if not x_set: raise ValueError("No compatible type!") from sage.libs.mpmath.all import call as mpcall from sage.libs.mpmath.all import chebyt as mpchebyt return mpcall(mpchebyt,args[0],args[-1],prec = precision) return mpcall(mpchebyt,args[0],x,parent=step_parent) def _maxima_init_evaled_(self, *args): """ EXAMPLES:: sage: chebyshev_T._maxima_init_evaled_(1,x) x 'x' sage: var('n') n sage: maxima(chebyshev_T(n,x)) chebyshev_t(n,x) """ n = args[0] x = args[1] return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x}) return maxima.eval('chebyshev_t({0},{1})'.format(n,x)) def _clenshaw_method_(self,*args): def _apply_formula_(self,*args): """ Clenshaw method for :class:`chebyshev_T` (means use recursions in this case). This is much faster for numerical evaluation than maxima! EXAMPLES:: sage: chebyshev_T._clenshaw_method_(2,0.1) == chebyshev_T._evalf_(2,0.1) sage: chebyshev_T._apply_formula_(2,0.1) == chebyshev_T._evalf_(2,0.1) True sage: chebyshev_T(51,x) 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 sage: chebyshev_T._clenshaw_method_(10,x) sage: chebyshev_T._apply_formula_(10,x) 512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1 """ help1 = 1 help2 = x if is_Expression(x) and k <= 50: if is_Expression(x) and k <= 25: # Recursion gives more compact representations for large k help3 = 0 for j in xrange(0,floor(k/2)+1): conversions=dict(maxima='chebyshev_u', mathematica='ChebyshevU')) def _clenshaw_method_(self,*args): def _apply_formula_(self,*args): """ Clenshaw method for :class:`chebyshev_U` (means we use the recursion) This is much faster for numerical evaluation than maxima. EXAMPLES:: sage: chebyshev_U._clenshaw_method_(2,0.1) == chebyshev_U._evalf_(2,0.1) sage: chebyshev_U._apply_formula_(2,0.1) == chebyshev_U._evalf_(2,0.1) True """ k = args[0] help1 = 1 help2 = 2*x if is_Expression(x) and k <= 50: if is_Expression(x) and k <= 25: # Recursion gives more compact representations for large k help3 = 0 for j in xrange(0,floor(k/2)+1): -2*((2*x + 1)*(2*x - 1)*x - 4*(2*x^2 - 1)*x)*(2*x + 1)*(2*x - 1) sage: abs(pari('polchebyshev(5, 2, 0.1)') - chebyshev_U(5,0.1)) < 1e-10 True """ if n == 0: EXAMPLES:: sage: x = PolynomialRing(QQ, 'x').gen() sage: chebyshev_T(2,x) 2*x^2 - 1 sage: maxima(chebyshev_U(5,x)) 32*x^5-32*x^3+6*x sage: var('n') n sage: maxima(chebyshev_U(n,x)) chebyshev_u(n,x) sage: maxima(chebyshev_U(2,x)) 4*x^2-1 """ n = args[0] x = args[1] return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x}) return maxima.eval('chebyshev_u({0},{1})'.format(n,x)) def _evalf_(self, *args,**kwds): """ Evaluate :class:`chebyshev_U` numerically with mpmath. If index is an integer use recursive formula since it is faster, for chebyshev polynomials. EXAMPLES:: 98280.0000000000 - 11310.0000000000*I sage: chebyshev_U(10,3).n(75) 4.661117900000000000000e7 sage: chebyshev_U._evalf_(1.5, Mod(8,9)) Traceback (most recent call last): ... ValueError: No compatible type! """ if args[0] in ZZ: return self._cheb_recur_(*args)[0] try: step_parent = kwds['parent'] real_parent = kwds['parent'] except KeyError: step_parent = parent(args[-1]) real_parent = parent(args[-1]) try: precision = step_parent.prec() except AttributeError: precision = RR.prec() x_set = False if hasattr(real_parent,"precision"): # Check if we have a data type with precision x = args[-1] step_parent = real_parent x_set = True else: if args[-1] in RR: x = RR(args[-1]) step_parent = RR x_set = True elif args[-1] in CC: x = CC(args[-1]) step_parent = CC x_set = True if not x_set: raise ValueError("No compatible type!") from sage.libs.mpmath.all import call as mpcall from sage.libs.mpmath.all import chebyu as mpchebyu return mpcall(mpchebyu,args[0],args[-1],prec = precision) return mpcall(mpchebyu,args[0],args[-1],parent = step_parent) def _eval_special_values_(self,*args): """ n + 1 sage: chebyshev_U(n,-1) (-1)^n*(n + 1) sage: chebyshev_U._eval_special_values_(26, Mod(0,9)) Traceback (most recent call last): ... ValueError: Value not found! sage: parent(chebyshev_U(3, Mod(8,9))) Ring of integers modulo 9 sage: parent(chebyshev_U(3, Mod(1,9))) Ring of integers modulo 9 sage: chebyshev_U(n, 0) 1/2*(-1)^(1/2*n)*((-1)^n + 1) sage: chebyshev_U(-3,x) + chebyshev_U(1,x) 0 sage: chebyshev_U._eval_special_values_(1.5, Mod(8,9)) Traceback (most recent call last): ... ValueError: No special values for non integral indices! """ if (not is_Expression(args[0])) and (not args[0] in ZZ): raise ValueError("No special values for non integral indices!") if args[-1] == 1: return (args[0]+1) return args[-1]*(args[0]+1) if args[-1] == -1: return (-1)**args[0]*(args[0]+1) return args[-1]**args[0]*(args[0]+1) if (args[-1] == 0): if (args[-1] == 0 and args[-1] in CC): return (1+(-1)**args[0])*(-1)**(args[0]/2)/2 if args[0] < 0 and args[0] in ZZ: return -self._eval_(-args[0]-2,args[-1]) raise ValueError("value not found") raise ValueError("Value not found!") def _eval_numpy_(self, *args): """