# 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


448  448  """ 
449  449  raise ValueError("no special values known") 
450  450  
 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^3x1) 
 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  
451  471  def _eval_(self, *args): 
452  472  """ 
453  473  The symbolic evaluation is either done with maxima, or with direct 
… 
… 

475  495  if not is_Expression(args[1]) and is_inexact(args[1]): 
476  496  try: 
477  497  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)) 
479  500  except AttributeError: 
480  501  pass 
481  502  except mpmath.NoConvergence: 
… 
… 

558  579  OrthogonalPolynomial.__init__(self, "chebyshev_T", nargs=2, 
559  580  conversions=dict(maxima='chebyshev_t', 
560  581  mathematica='ChebyshevT')) 
561   
 582  
562  583  def _eval_special_values_(self,*args): 
563  584  """ 
564  585  Values known for special values of x. 
… 
… 

596  617  3363.00000000000 
597  618  sage: chebyshev_T(5,0.3).n() 
598  619  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 
599  622  """ 
600  623  try: 
601  624  step_parent = kwds['parent'] 
… 
… 

625  648  x = args[1] 
626  649  return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x}) 
627  650  
 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=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  
628  673  def _clenshaw_method_(self,*args): 
629  674  """ 
630  675  Clenshaw method for :class:`chebyshev_T` (means use recursions in this 
… 
… 

636  681  sage: chebyshev_T._clenshaw_method_(2,0.1) == chebyshev_T._evalf_(2,0.1) 
637  682  True 
638  683  """ 
639   k = args[0] 
 684  n = args[0] 
640  685  x = args[1] 
641  686  
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 
646  692  
647  693  # TODO: When evaluation of Symbolic Expressions works better 
648  694  # use these explicit formulas instead! 
… 
… 

652  698  # return cosh(k*acosh(x)) 
653  699  #else: # x < 1 
654  700  # 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(kj1) / factorial(j) / factorial(k2*j) 
663   help3 = help3 + (1)**j * (2*x)**(k2*j) * f 
664   help3 = help3 * k / 2 
665   else: 
666   for j in xrange(0,k1): 
667   help3 = 2*x*help2  help1 
668   help1 = help2 
669   help2 = help3 
670   
671   return help3 
 701  
672  702  
673  703  def _eval_numpy_(self, *args): 
674  704  """ 
… 
… 

709  739  sage: derivative(chebyshev_T(k,x),k) 
710  740  Traceback (most recent call last): 
711  741  ... 
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 
713  743  """ 
714  744  diff_param = kwds['diff_param'] 
715  745  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") 
717  747  
718  748  return args[0]*chebyshev_U(args[0]1,args[1]) 
719  749  
… 
… 

749  779  OrthogonalPolynomial.__init__(self, "chebyshev_U", nargs=2, 
750  780  conversions=dict(maxima='chebyshev_u', 
751  781  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=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**21 
 800  
 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] 
752  807  
753  808  def _clenshaw_method_(self,*args): 
754  809  """ 
… 
… 

763  818  """ 
764  819  k = args[0] 
765  820  x = args[1] 
766   
767   if k == 0: 
768   return 1 
769   if k == 1: 
770   return 2*x 
771  821  
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(kj) / factorial(j) / factorial(k2*j) # Change to a binomial? 
779   help3 = help3 + (1)**j * (2*x)**(k2*j) * f 
780   
781   else: 
782   for j in xrange(0,k1): 
783   help3 = 2*x*help2  help1 
784   help1 = help2 
785   help2 = help3 
786   
787   return help3 
 822  return self._cheb_recur_(k,x)[0] 
788  823  
789  824  def _maxima_init_evaled_(self, *args): 
790  825  """ 
… 
… 

818  853  
819  854  try: 
820  855  precision = step_parent.prec() 
821   except AttributeError: 
 856  except ValueError: 
822  857  precision = RR.prec() 
823  858  
824  859  from sage.libs.mpmath.all import call as mpcall 
… 
… 

890  925  sage: derivative(chebyshev_U(k,x),k) 
891  926  Traceback (most recent call last): 
892  927  ... 
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 
894  929  """ 
895  930  diff_param = kwds['diff_param'] 
896  931  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") 
898  933  
899  934  return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])args[1]* 
900  935  chebyshev_U(args[0],args[1]))/(args[1]**21) 