# 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


315  315  # http://www.gnu.org/licenses/ 
316  316  #***************************************************************************** 
317  317  
 318  import warnings 
 319  
318  320  from sage.misc.sage_eval import sage_eval 
319   from sage.rings.all import ZZ, RR 
 321  from sage.rings.all import ZZ, RR, CC, RIF, CIF 
320  322  from sage.calculus.calculus import maxima 
321  323  
322  324  
… 
… 

399  401  except KeyError: 
400  402  self._maxima_name = None 
401  403  
402   BuiltinFunction.__init__(self, name=name, nargs=nargs, 
 404  super(OrthogonalPolynomial,self).__init__(name=name, nargs=nargs, 
403  405  latex_name=latex_name, conversions=conversions) 
404  406  
405  407  def _maxima_init_evaled_(self, *args): 
… 
… 

420  422  """ 
421  423  return None 
422  424  
423   def _clenshaw_method_(self, *args): 
 425  def _apply_formula_(self, *args): 
424  426  """ 
425  427  The Clenshaw method uses the three term recursion of the polynomial, 
426  428  or explicit formulas instead of maxima to evaluate the polynomial 
… 
… 

432  434  
433  435  sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial 
434  436  sage: new = OrthogonalPolynomial('testo_P') 
435   sage: new._clenshaw_method_(1,2.0) 
 437  sage: new._apply_formula_(1,2.0) 
436  438  Traceback (most recent call last): 
437  439  ... 
438  440  NotImplementedError: no recursive calculation of values implemented 
… 
… 

474  476  64*x^7  112*x^5 + 56*x^3  7*x 
475  477  sage: chebyshev_T(3/2,x) 
476  478  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 
477  484  sage: parent(chebyshev_T(4, RIF(5))) 
478  485  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 
479  496  """ 
480  497  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): 
482  500  try: 
483  501  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) 
486  503  except AttributeError: 
487  504  pass 
488  505  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) 
491  508  except ValueError: 
492  509  pass 
493  510  
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]): 
502  514  return None 
 515  
 516  # Check for known identities 
503  517  try: 
504  518  return self._eval_special_values_(*args) 
505  519  except ValueError: 
… 
… 

511  525  if args[0] < 0 and args[0] in ZZ: 
512  526  return None 
513  527  
514   if not is_Expression(args[0]): 
 528  if args[0] in ZZ: 
515  529  try: 
516   return self._clenshaw_method_(*args) 
 530  return self._apply_formula_(*args) 
517  531  except NotImplementedError: 
518  532  pass 
519  533  
520  534  if self._maxima_name is None: 
521  535  return None 
522  536  
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: 
528  540  return None 
529  541  
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): 
536  543  """ 
537  544  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 
539  546  would normally allow. 
540  547  Thus we have the distinction between algebraic objects (if n is an integer), 
541  548  and else as symbolic function. 
… 
… 

552  559  sage: parent(chebyshev_T(5, x)) 
553  560  Univariate Polynomial Ring in x over Rational Field 
554  561  """ 
 562  if 'hold' not in kwds: 
 563  kwds['hold'] = False 
 564  if 'coerce' not in kwds: 
 565  kwds['coerce']=True 
555  566  
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 
558  569  else: # Consider OrthogonalPolynomial as symbol 
559   return super(OrthogonalPolynomial,self).__call__(n,x) 
 570  return super(OrthogonalPolynomial,self).__call__(*args,**kwds) 
560  571  
 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 
561  584  
562  585  class Func_chebyshev_T(OrthogonalPolynomial): 
563  586  """ 
… 
… 

590  613  sage: chebyshev_T2(1,x) 
591  614  x 
592  615  """ 
593   OrthogonalPolynomial.__init__(self, "chebyshev_T", nargs=2, 
 616  super(Func_chebyshev_T,self).__init__("chebyshev_T", nargs=2, 
594  617  conversions=dict(maxima='chebyshev_t', 
595  618  mathematica='ChebyshevT')) 
596  619  
… 
… 

607  630  1 
608  631  sage: chebyshev_T(n,1) 
609  632  (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! 
610  647  """ 
 648  if (not is_Expression(args[0])) and (not args[0] in ZZ): 
 649  raise ValueError("No special values for non integral indices!") 
 650  
611  651  if args[1] == 1: 
612   return 1 
 652  return args[1] 
613  653  
614  654  if args[1] == 1: 
615   return (1)**args[0] 
 655  return args[1]**args[0] 
616  656  
617   if (args[1] == 0): 
 657  if (args[1] == 0 and args[1] in CC): 
618  658  return (1+(1)**args[0])*(1)**(args[0]/2)/2 
619  659  
620  660  if args[0] < 0 and args[0] in ZZ: 
621  661  return self._eval_(args[0],args[1]) 
622  662  
623   raise ValueError("value not found") 
 663  raise ValueError("Value not found!") 
624  664  
625  665  def _evalf_(self, *args,**kwds): 
626  666  """ 
627  667  Evaluates :class:`chebyshev_T` numerically with mpmath. 
 668  If the index is an integer we use the recursive formula since 
 669  it is faster. 
628  670  
629  671  EXAMPLES:: 
630  672  
… 
… 

636  678  0.998880000000000 
637  679  sage: chebyshev_T(1/2, 0) 
638  680  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  
639  686  """ 
 687  if args[0] in ZZ: 
 688  return self._cheb_recur_(*args)[0] 
 689  
640  690  try: 
641   step_parent = kwds['parent'] 
 691  real_parent = kwds['parent'] 
642  692  except KeyError: 
643   step_parent = parent(args[1]) 
 693  real_parent = parent(args[1]) 
644  694  
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 
649  709  
 710  if not x_set: 
 711  raise ValueError("No compatible type!") 
 712  
650  713  from sage.libs.mpmath.all import call as mpcall 
651  714  from sage.libs.mpmath.all import chebyt as mpchebyt 
652  715  
653   return mpcall(mpchebyt,args[0],args[1],prec = precision) 
 716  return mpcall(mpchebyt,args[0],x,parent=step_parent) 
654  717  
655  718  def _maxima_init_evaled_(self, *args): 
656  719  """ 
… 
… 

659  722  EXAMPLES:: 
660  723  
661  724  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  
663  731  """ 
664  732  n = args[0] 
665  733  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)) 
667  735  
668   def _clenshaw_method_(self,*args): 
 736  
 737  def _apply_formula_(self,*args): 
669  738  """ 
670  739  Clenshaw method for :class:`chebyshev_T` (means use recursions in this 
671  740  case). This is much faster for numerical evaluation than maxima! 
… 
… 

674  743  
675  744  EXAMPLES:: 
676  745  
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) 
678  747  True 
679  748  sage: chebyshev_T(51,x) 
680  749  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) 
682  751  512*x^10  1280*x^8 + 1120*x^6  400*x^4 + 50*x^2  1 
683  752  
684  753  """ 
… 
… 

701  770  
702  771  help1 = 1 
703  772  help2 = x 
704   if is_Expression(x) and k <= 50: 
 773  if is_Expression(x) and k <= 25: 
705  774  # Recursion gives more compact representations for large k 
706  775  help3 = 0 
707  776  for j in xrange(0,floor(k/2)+1): 
… 
… 

814  883  conversions=dict(maxima='chebyshev_u', 
815  884  mathematica='ChebyshevU')) 
816  885  
817   def _clenshaw_method_(self,*args): 
 886  def _apply_formula_(self,*args): 
818  887  """ 
819  888  Clenshaw method for :class:`chebyshev_U` (means we use the recursion) 
820  889  This is much faster for numerical evaluation than maxima. 
… 
… 

823  892  
824  893  EXAMPLES:: 
825  894  
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) 
827  896  True 
828  897  """ 
829  898  k = args[0] 
… 
… 

836  905  
837  906  help1 = 1 
838  907  help2 = 2*x 
839   if is_Expression(x) and k <= 50: 
 908  if is_Expression(x) and k <= 25: 
840  909  # Recursion gives more compact representations for large k 
841  910  help3 = 0 
842  911  for j in xrange(0,floor(k/2)+1): 
… 
… 

862  931  2*((2*x + 1)*(2*x  1)*x  4*(2*x^2  1)*x)*(2*x + 1)*(2*x  1) 
863  932  sage: abs(pari('polchebyshev(5, 2, 0.1)')  chebyshev_U(5,0.1)) < 1e10 
864  933  True 
865   
866   
867  934  """ 
868  935  
869  936  if n == 0: 
… 
… 

883  950  
884  951  EXAMPLES:: 
885  952  
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^532*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^21 
889  961  """ 
890  962  n = args[0] 
891  963  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)) 
893  965  
894  966  def _evalf_(self, *args,**kwds): 
895  967  """ 
896  968  Evaluate :class:`chebyshev_U` numerically with mpmath. 
 969  If index is an integer use recursive formula since it is faster, 
 970  for chebyshev polynomials. 
897  971  
898  972  EXAMPLES:: 
899  973  
… 
… 

901  975  98280.0000000000  11310.0000000000*I 
902  976  sage: chebyshev_U(10,3).n(75) 
903  977  4.661117900000000000000e7 
 978  sage: chebyshev_U._evalf_(1.5, Mod(8,9)) 
 979  Traceback (most recent call last): 
 980  ... 
 981  ValueError: No compatible type! 
904  982  """ 
 983  if args[0] in ZZ: 
 984  return self._cheb_recur_(*args)[0] 
905  985  try: 
906   step_parent = kwds['parent'] 
 986  real_parent = kwds['parent'] 
907  987  except KeyError: 
908   step_parent = parent(args[1]) 
 988  real_parent = parent(args[1]) 
909  989  
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!") 
914  1007  
915  1008  from sage.libs.mpmath.all import call as mpcall 
916  1009  from sage.libs.mpmath.all import chebyu as mpchebyu 
917  1010  
918   return mpcall(mpchebyu,args[0],args[1],prec = precision) 
 1011  return mpcall(mpchebyu,args[0],args[1],parent = step_parent) 
919  1012  
920  1013  def _eval_special_values_(self,*args): 
921  1014  """ 
… 
… 

929  1022  n + 1 
930  1023  sage: chebyshev_U(n,1) 
931  1024  (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! 
932  1041  """ 
 1042  if (not is_Expression(args[0])) and (not args[0] in ZZ): 
 1043  raise ValueError("No special values for non integral indices!") 
 1044  
933  1045  if args[1] == 1: 
934   return (args[0]+1) 
 1046  return args[1]*(args[0]+1) 
935  1047  
936  1048  if args[1] == 1: 
937   return (1)**args[0]*(args[0]+1) 
 1049  return args[1]**args[0]*(args[0]+1) 
938  1050  
939   if (args[1] == 0): 
 1051  if (args[1] == 0 and args[1] in CC): 
940  1052  return (1+(1)**args[0])*(1)**(args[0]/2)/2 
941  1053  
942  1054  if args[0] < 0 and args[0] in ZZ: 
943  1055  return self._eval_(args[0]2,args[1]) 
944  1056  
945   raise ValueError("value not found") 
 1057  raise ValueError("Value not found!") 
946  1058  
947  1059  def _eval_numpy_(self, *args): 
948  1060  """ 