Ticket #9706: orthogonal_polys.8.py

File orthogonal_polys.8.py, 80.6 KB (added by maldun, 11 years ago)

Latest version with some code cleanup (no program changes)

Line 
1# -*- coding: utf-8 -*-
2r"""
3Orthogonal Polynomials
4
5This module wraps some of the orthogonal/special functions in the
6Maxima package "orthopoly". This package was written by Barton
7Willis of the University of Nebraska at Kearney. It is released
8under the terms of the General Public License (GPL). Send
9Maxima-related bug reports and comments on this module to
10willisb@unk.edu. In your report, please include Maxima and specfun
11version information.
12-----------------------------------------------------------------
13Update (2010): The orthogonal polynomials are now selfstanding classes.
14Altough some of the use the old maxima wrappers, most of the evaluation
15is done with the help of direct methods. Also the orthogonal polys
16are handled as symbolic expressions, if the input is not known
17
18
19-  The Chebyshev polynomial of the first kind arises as a solution
20   to the differential equation
21   
22   .. math::
23
24         (1-x^2)\,y'' - x\,y' + n^2\,y = 0
25
26
27   and those of the second kind as a solution to
28   
29   .. math::
30
31         (1-x^2)\,y'' - 3x\,y' + n(n+2)\,y = 0.
32
33
34   The Chebyshev polynomials of the first kind are defined by the
35   recurrence relation
36   
37   .. math::
38
39     T_0(x) = 1 \, T_1(x) = x \, T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x). \,
40
41
42   The Chebyshev polynomials of the second kind are defined by the
43   recurrence relation
44   
45   .. math::
46
47     U_0(x) = 1 \, U_1(x) = 2x \, U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x). \,
48
49
50
51   For integers `m,n`, they satisfy the orthogonality
52   relations
53   
54   .. math::
55
56     \int_{-1}^1 T_n(x)T_m(x)\,\frac{dx}{\sqrt{1-x^2}} =\left\{ \begin{matrix} 0 &: n\ne m~~~~~\\ \pi &: n=m=0\\ \pi/2 &: n=m\ne 0 \end{matrix} \right. 
57
58
59   and
60
61   
62   .. math::
63
64     \int_{-1}^1 U_n(x)U_m(x)\sqrt{1-x^2}\,dx =\frac{\pi}{2}\delta_{m,n}.
65
66
67
68   They are named after Pafnuty Chebyshev (alternative
69   transliterations: Tchebyshef or Tschebyscheff).
70
71-  The Hermite polynomials are defined either by
72   
73   .. math::
74
75     H_n(x)=(-1)^n e^{x^2/2}\frac{d^n}{dx^n}e^{-x^2/2}
76
77
78   (the "probabilists' Hermite polynomials"), or by
79
80   
81   .. math::
82
83     H_n(x)=(-1)^n e^{x^2}\frac{d^n}{dx^n}e^{-x^2}
84
85
86   (the "physicists' Hermite polynomials"). Sage (via Maxima)
87   implements the latter flavor. These satisfy the orthogonality
88   relation
89   
90   .. math::
91
92     \int_{-\infty}^\infty H_n(x)H_m(x)\,e^{-x^2}\,dx ={n!2^n}{\sqrt{\pi}}\delta_{nm}
93
94
95
96   They are named in honor of Charles Hermite.
97
98-  Each *Legendre polynomial* `P_n(x)` is an `n`-th degree polynomial.
99   It may be expressed using Rodrigues' formula:
100
101   .. math::
102
103      P_n(x) = (2^n n!)^{-1} {\frac{d^n}{dx^n} } \left[ (x^2 -1)^n \right].
104
105   These are solutions to Legendre's differential equation:
106
107   .. math::
108
109      {\frac{d}{dx}} \left[ (1-x^2) {\frac{d}{dx}} P(x) \right] + n(n+1)P(x) = 0.
110
111   and satisfy the orthogonality relation
112
113   .. math::
114   
115      \int_{-1}^{1} P_m(x) P_n(x)\,dx = {\frac{2}{2n + 1}} \delta_{mn}
116
117   The *Legendre function of the second kind* `Q_n(x)` is another
118   (linearly independent) solution to the Legendre differential equation.
119   It is not an "orthogonal polynomial" however.
120
121   The associated Legendre functions of the first kind
122   `P_\ell^m(x)` can be given in terms of the "usual"
123   Legendre polynomials by
124
125   .. math::
126
127     \begin{array}{ll} P_\ell^m(x)    &=  (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}P_\ell(x) \\ &=  \frac{(-1)^m}{2^\ell \ell!} (1-x^2)^{m/2}\frac{d^{\ell+m}}{dx^{\ell+m}}(x^2-1)^\ell. \end{array}
128
129
130   Assuming `0 \le m \le \ell`, they satisfy the orthogonality
131   relation:
132   
133   .. math::
134
135      \int_{-1}^{1} P_k ^{(m)} P_\ell ^{(m)} dx  = \frac{2 (\ell+m)!}{(2\ell+1)(\ell-m)!}\ \delta _{k,\ell},
136
137
138   where `\delta _{k,\ell}` is the Kronecker delta.
139
140   The associated Legendre functions of the second kind
141   `Q_\ell^m(x)` can be given in terms of the "usual"
142   Legendre polynomials by
143
144   
145   .. math::
146
147     Q_\ell^m(x)   =  (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}Q_\ell(x).
148
149
150
151   They are named after Adrien-Marie Legendre.
152
153-  Laguerre polynomials may be defined by the Rodrigues formula
154   
155   .. math::
156
157      L_n(x)=\frac{e^x}{n!}\frac{d^n}{dx^n}\left(e^{-x} x^n\right).
158
159
160   They are solutions of Laguerre's equation:
161
162   
163   .. math::
164
165      x\,y'' + (1 - x)\,y' + n\,y = 0\, 
166
167   and satisfy the orthogonality relation
168
169   
170   .. math::
171
172      \int_0^\infty L_m(x) L_n(x) e^{-x}\,dx = \delta_{mn}.
173
174
175
176   The generalized Laguerre polynomials may be defined by the
177   Rodrigues formula:
178
179   
180   .. math::
181
182       L_n^{(\alpha)}(x)   = {\frac{x^{-\alpha} e^x}{n!}}{\frac{d^n}{dx^n}} \left(e^{-x} x^{n+\alpha}\right) .
183
184
185   (These are also sometimes called the associated Laguerre
186   polynomials.) The simple Laguerre polynomials are recovered from
187   the generalized polynomials by setting `\alpha =0`.
188
189   They are named after Edmond Laguerre.
190
191-  Jacobi polynomials are a class of orthogonal polynomials. They
192   are obtained from hypergeometric series in cases where the series
193   is in fact finite:
194   
195   .. math::
196
197     P_n^{(\alpha,\beta)}(z) =\frac{(\alpha+1)_n}{n!} \,_2F_1\left(-n,1+\alpha+\beta+n;\alpha+1;\frac{1-z}{2}\right) ,
198
199
200   where `()_n` is Pochhammer's symbol (for the rising
201   factorial), (Abramowitz and Stegun p561.) and thus have the
202   explicit expression
203
204   
205   .. math::
206
207     P_n^{(\alpha,\beta)} (z) = \frac{\Gamma (\alpha+n+1)}{n!\Gamma (\alpha+\beta+n+1)} \sum_{m=0}^n {n\choose m} \frac{\Gamma (\alpha + \beta + n + m + 1)}{\Gamma (\alpha + m + 1)} \left(\frac{z-1}{2}\right)^m . 
208
209
210
211   They are named after Carl Jacobi.
212
213-  Ultraspherical or Gegenbauer polynomials are given in terms of
214   the Jacobi polynomials `P_n^{(\alpha,\beta)}(x)` with
215   `\alpha=\beta=a-1/2` by
216
217   
218   .. math::
219
220     C_n^{(a)}(x)= \frac{\Gamma(a+1/2)}{\Gamma(2a)}\frac{\Gamma(n+2a)}{\Gamma(n+a+1/2)} P_n^{(a-1/2,a-1/2)}(x).
221
222
223   They satisfy the orthogonality relation
224   
225   .. math::
226
227     \int_{-1}^1(1-x^2)^{a-1/2}C_m^{(a)}(x)C_n^{(a)}(x)\, dx =\delta_{mn}2^{1-2a}\pi \frac{\Gamma(n+2a)}{(n+a)\Gamma^2(a)\Gamma(n+1)} ,
228
229
230   for `a>-1/2`. They are obtained from hypergeometric series
231   in cases where the series is in fact finite:
232
233   
234   .. math::
235
236     C_n^{(a)}(z) =\frac{(2a)^{\underline{n}}}{n!} \,_2F_1\left(-n,2a+n;a+\frac{1}{2};\frac{1-z}{2}\right)
237
238
239   where `\underline{n}` is the falling factorial. (See
240   Abramowitz and Stegun p561)
241
242   They are named for Leopold Gegenbauer (1849-1903).
243
244
245For completeness, the Pochhammer symbol, introduced by Leo August
246Pochhammer, `(x)_n`, is used in the theory of special
247functions to represent the "rising factorial" or "upper factorial"
248
249.. math::
250
251         (x)_n=x(x+1)(x+2)\cdots(x+n-1)=\frac{(x+n-1)!}{(x-1)!}.
252
253
254On the other hand, the "falling factorial" or "lower factorial" is
255
256.. math::
257
258     x^{\underline{n}}=\frac{x!}{(x-n)!} ,
259
260
261in the notation of Ronald L. Graham, Donald E. Knuth and Oren
262Patashnik in their book Concrete Mathematics.
263
264.. note::
265
266   The first call of any of these will usually cost a bit extra
267   (it loads "specfun", but I'm not sure if that is the real reason).
268   The next call is usually faster but not always.
269
270TODO: Implement associated Legendre polynomials and Zernike
271polynomials. (Neither is in Maxima.)
272http://en.wikipedia.org/wiki/Associated_Legendre_polynomials
273http://en.wikipedia.org/wiki/Zernike_polynomials
274
275REFERENCES:
276
277-  Abramowitz and Stegun: Handbook of Mathematical Functions,
278   http://www.math.sfu.ca/ cbm/aands/
279
280-  http://en.wikipedia.org/wiki/Chebyshev_polynomials
281
282-  http://en.wikipedia.org/wiki/Legendre_polynomials
283
284-  http://en.wikipedia.org/wiki/Hermite_polynomials
285
286-  http://mathworld.wolfram.com/GegenbauerPolynomial.html
287
288-  http://en.wikipedia.org/wiki/Jacobi_polynomials
289
290-  http://en.wikipedia.org/wiki/Laguerre_polynomia
291
292-  http://en.wikipedia.org/wiki/Associated_Legendre_polynomials
293
294
295AUTHORS:
296
297- David Joyner (2006-06)
298- Stefan Reiterer (2010-)
299"""
300
301#*****************************************************************************
302#       Copyright (C) 2006 William Stein <wstein@gmail.com>
303#                     2006 David Joyner <wdj@usna.edu>
304#                     2010 Stefan Reiterer <stefan.reiterer@uni-graz.at>
305#
306#  Distributed under the terms of the GNU General Public License (GPL)
307#
308#    This code is distributed in the hope that it will be useful,
309#    but WITHOUT ANY WARRANTY; without even the implied warranty of
310#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
311#    General Public License for more details.
312#
313#  The full text of the GPL is available at:
314#
315#                  http://www.gnu.org/licenses/
316#*****************************************************************************
317
318import copy
319import sage.plot.plot
320import sage.interfaces.all
321from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
322from sage.rings.rational_field import RationalField
323from sage.rings.real_mpfr import RealField
324from sage.misc.sage_eval import sage_eval
325from sage.rings.all import QQ, ZZ, RR, CDF, RDF
326from sage.symbolic.ring import SR
327import sage.rings.commutative_ring as commutative_ring
328import sage.rings.ring as ring
329from sage.calculus.calculus import maxima
330from sage.symbolic.function import BuiltinFunction, GinacFunction, is_inexact
331from sage.symbolic.expression import is_Expression
332import sage.functions.special
333from sage.functions.special import MaximaFunction, meval
334import numpy,scipy
335import math
336#import sage.libs.mpmath.all as mpmath #This line does evil things...
337from sage.functions.other import floor, gamma, factorial, abs, binomial
338from sage.functions.other import sqrt, conjugate
339from sage.functions.trig import sin, cos
340from sage.functions.log import ln
341import sage.symbolic.expression as expression
342from sage.structure.parent import Parent
343from sage.structure.coerce import parent
344
345_done = False
346def _init():
347    """
348    Internal function which checks if Maxima has loaded the
349    "orthopoly" package.  All functions using this in this
350    file should call this function first.
351   
352    TEST:
353
354    The global starts ``False``::
355
356        sage: sage.functions.orthogonal_polys._done
357        False
358
359    Then after using one of these functions, it changes:: (The value is now
360    False for chebyshev_T because chebyshev_T uses clenshaw method instead...)
361
362        sage: from sage.functions.orthogonal_polys import chebyshev_T
363        sage: chebyshev_T(2,x)
364        2*x^2 - 1
365        sage: sage.functions.orthogonal_polys._done
366        False
367        sage: legendre_Q(1,x)
368        1/2*x*log(-(x + 1)/(x - 1)) - 1
369        sage: sage.functions.orthogonal_polys._done
370        True
371   
372    Note that because here we use a Pynac variable ``x``,
373    the representation of the function is different from
374    its actual doctest, where a polynomial indeterminate
375    ``x`` is used.
376    """
377    global _done
378    if _done:
379        return
380    maxima.eval('load("orthopoly");')
381    # TODO -- make it possible to use the intervals returned
382    # instead of just discarding this info!
383    maxima.eval('orthopoly_returns_intervals:false;')
384    _done = True
385
386#load /home/maldun/sage/sage-4.5.2/devel/sage-ortho/sage/functions/orthogonal_polys.py
387
388class OrthogonalPolynomial(BuiltinFunction):
389    """
390    Base Class for Orthogonal Polynomials. The evaluation as a polynomial
391    is done via maxima due performance reasons. Therefore the internal name
392    in maxima maxima_name has to be declared.
393    Convention: The first argument is always the order of the polynomial,
394    he last one is always the value x where the polynomial is evaluated.
395
396    """
397    def __init__(self, name, nargs = 2, latex_name = None, conversions = {}):
398        try:
399            self._maxima_name = conversions['maxima'] 
400        except KeyError:
401            self._maxima_name = None
402
403        BuiltinFunction.__init__(self, name = name,
404        nargs = nargs, latex_name = latex_name, conversions = conversions)
405
406    def _maxima_init_evaled_(self, *args):
407        """
408        Returns a string which represents this function evaluated at
409        *args* in Maxima.
410        In fact these are thought to be the old wrappers for the orthogonal
411        polynomials. These are used when the other evaluation methods fail,
412        or are not fast enough. Experiments showed that for the symbolic
413        evaluation for larger n maxima is faster, but for small n simply use
414        of the recursion formulas is faster. A little switch does the trick...
415
416        EXAMPLES::
417
418            sage: chebyshev_T(3,x)
419            4*x^3 - 3*x
420        """
421        return None
422
423    def _clenshaw_method_(self,*args):
424        """
425        The Clenshaw method uses the three term recursion of the polynomial,
426        or explicit formulas instead of maxima to evaluate the polynomial
427        efficiently, if the x argument is not a symbolic expression.
428        The name comes from the Clenshaw algorithm for fast evaluation of
429        polynomialsin chebyshev form.
430       
431        comparison: Maxima and Clenshaw algorithm for non-symbolic evaluation:
432         #sage: time chebyshev_T(50,10) #clenshaw
433         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
434         #Wall time: 0.00 s
435         #49656746733678312490954442369580252421769338391329426325400124999
436         #sage: time sage.functions.orthogonal_polys.chebyshev_T(50,10)
437         #maxima
438         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
439         #Wall time: 0.05 s
440         #49656746733678312490954442369580252421769338391329426325400124999
441
442         #sage: time chebyshev_T(500,10); #clenshaw
443         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
444         #Wall time: 0.00 s
445         #sage: time sage.functions.orthogonal_polys.chebyshev_T(500,10);
446         #maxima
447         #CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s
448         #Wall time: 0.77 s
449        """
450        raise NotImplementedError(
451            "No recursive calculation of values implemented (yet)!")
452
453    def _eval_special_values_(self,*args):
454        """
455        Evals the polynomial explicitly for special values.
456        EXAMPLES:
457           
458            sage: var('n')
459            n
460            sage: chebyshev_T(n,-1)
461            (-1)^n
462        """
463        raise ValueError("No special values known!")
464           
465
466    def _eval_(self, *args):
467        """
468       
469        The symbolic evaluation is done with maxima, because the evaluation of
470        the Polynomial representation seems to be quite clever.
471        For the fast numerical evaluation an other method should be used...
472        Therefore I suggest Clenshaw's algorithm, which uses the rekursion!
473        The function also checks for special values, and if
474        the order is an integer and in range!
475       
476        performance:
477        #sage: time chebyshev_T(5,x) #maxima
478        #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
479        #Wall time: 0.16 s
480        #16*(x - 1)^5 + 80*(x - 1)^4 + 140*(x - 1)^3 + 100*(x - 1)^2 + 25*x - 24
481
482        #sage: time chebyshev_T(5,x) #clenshaw
483        #CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
484        #Wall time: 0.01 s
485        #16*x^5 - 20*x^3 + 5*x
486
487        #time chebyshev_T(50,x)
488        #CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
489        #Wall time: 0.04 s
490        #562949953421312*(x - 1)^50 + 28147497671065600*(x - 1)^49 +....
491
492        #time chebyshev_T(100,x);
493        #CPU times: user 0.08 s, sys: 0.00 s, total: 0.08 s
494        #Wall time: 0.08 s
495
496        EXAMPLES::
497            sage: chebyshev_T(5,x)
498            16*x^5 - 20*x^3 + 5*x 
499            sage: var('n')
500            n
501            sage: chebyshev_T(n,-1)
502            (-1)^n
503            sage: chebyshev_T(-7,x)
504            chebyshev_T(-7, x)
505            sage: chebyshev_T(3/2,x)
506            chebyshev_T(3/2, x)
507   
508        """
509       
510        if not is_Expression(args[0]):
511           
512            if not is_Expression(args[-1]) and is_inexact(args[-1]):
513                try:
514                    import sage.libs.mpmath.all as mpmath
515                    return self._evalf_(*args) 
516                except AttributeError:
517                    pass
518                except mpmath.NoConvergence:
519                    print "Warning: mpmath returns NoConvergence!" 
520                    print "Switching to clenshaw_method, but it \
521                    may not be stable!"
522                except ValueError:
523                    pass
524
525            #A faster check would be nice...
526            if args[0] != floor(args[0]):
527                if not is_Expression(args[-1]):
528                    try:
529                        return self._evalf_(*args)
530                    except AttributeError:
531                        pass
532                else:
533                    return None
534       
535            if args[0] < 0:
536                return None
537           
538           
539        try:
540            return self._eval_special_values_(*args)
541        except ValueError:
542            pass
543
544   
545        if not is_Expression(args[0]):
546           
547            try: 
548                return self._clenshaw_method_(*args)
549            except NotImplementedError:
550                pass
551       
552        if self._maxima_name is None:
553            return None
554        else:
555            _init()
556            try:
557                #s = maxima(self._maxima_init_evaled_(*args))
558                #This above is very inefficient! The older
559                #methods were much faster...
560                return self._maxima_init_evaled_(*args)
561            except TypeError:
562                return None
563            if self._maxima_name in repr(s):
564                return None
565            else:
566                return s.sage()
567
568    #def _eval_numpy_(self, *args):
569        #TODO: numpy_eval with help of the new scipy version!!!!
570        #Reason scipy 0.8 supports stable and fast numerical evaluation
571        #of ortho polys
572        #Now this only evaluates the array pointwise, and only the first one...
573        #if isinstance(arg[0], numpy.ndarray)
574       
575       
576       
577       
578class Func_chebyshev_T(OrthogonalPolynomial): 
579
580    """
581    Class for the Chebyshev polynomial of the first kind.
582   
583    REFERENCE:
584
585    - AS 22.5.31 page 778 and AS 6.1.22 page 256.
586   
587    EXAMPLES::
588   
589       sage: chebyshev_T(3,x)
590       4*x^3 - 3*x
591       sage: var('k')
592       k
593       sage: test = chebyshev_T(k,x)
594       sage: test
595       chebyshev_T(k, x)
596    """
597   
598    def __init__(self):
599        OrthogonalPolynomial.__init__(self,"chebyshev_T",nargs = 2,
600conversions =dict(maxima='chebyshev_t',mathematica='ChebyshevT'))
601   
602    def _eval_special_values_(self,*args):
603        """
604        Values known for special values of x.
605        For details see A.S. 22.4 (p. 777)
606       
607        EXAMPLES:
608           
609            sage: var('n')
610            n
611            sage: chebyshev_T(n,1)
612            1
613            sage: chebyshev_T(n,-1)
614            (-1)^n
615        """
616        if args[-1] == 1:
617            return 1
618       
619        if args[-1] == -1:
620            return (-1)**args[0]
621
622        if (args[-1] == 0):
623            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
624
625        raise ValueError("Value not found!")
626   
627    def _evalf_(self, *args,**kwds):
628        """
629        Evals chebyshev_T
630        numerically with mpmath.
631        EXAMPLES::
632            sage: chebyshev_T(10,3).n(75)
633            2.261953700000000000000e7
634        """
635        try:
636            step_parent = kwds['parent']
637        except KeyError:
638            step_parent = parent(args[-1])
639
640        import sage.libs.mpmath.all as mpmath
641
642        try:
643            precision = step_parent.prec()
644        except AttributeError:
645            precision = mpmath.mp.prec
646
647        return mpmath.call(mpmath.chebyt,args[0],args[-1],prec = precision)
648
649    def _maxima_init_evaled_(self, *args):
650        n = args[0]
651        x = args[1]
652        return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
653
654    def _clenshaw_method_(self,*args):
655            """
656            Clenshaw method for chebyshev_T (means use recursions in this case)
657            This is much faster for numerical evaluation than maxima!
658            See A.S. 227 (p. 782) for details for the recurions
659            """
660
661            k = args[0]
662            x = args[1]
663
664            if k == 0:
665                return 1
666            elif k == 1:
667                return x
668            else:
669                #TODO: When evaluation of Symbolic Expressions works better
670                #use these explicit formulas instead!
671                #if -1 <= x <= 1:
672                #    return cos(k*acos(x))
673                #elif 1 < x:
674                #    return cosh(k*acosh(x))
675                #else: # x < -1
676                #    return (-1)**(k%2)*cosh(k*acosh(-x))
677               
678                help1 = 1
679                help2 = x
680                if is_Expression(x):
681                    #raise NotImplementedError
682                    help3 = 0
683                    for j in xrange(0,floor(k/2)+1):
684                        help3 = \
685help3 +(-1)**j*(2*x)**(k-2*j)*factorial(k-j-1)/factorial(j)/factorial(k-2*j)
686                    help3 = help3*k/2
687                else:
688                     for j in xrange(0,k-1):
689                        help3 = 2*x*help2 - help1
690                        help1 = help2
691                        help2 = help3
692               
693                return help3
694
695    def _eval_numpy_(self, *args):
696        #TODO: numpy_eval with help of the new scipy version!!!!
697        #Reason scipy 0.8 supports stable and fast numerical evaluation
698        #of ortho polys
699        #Now this only evaluates the array pointwise, and only the first one...
700        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
701        """
702        EXAMPLES::
703        sage: import numpy
704        sage: z = numpy.array([1,2])
705        sage: z2 = numpy.array([[1,2],[1,2]])
706        sage: z3 = numpy.array([1,2,3.])
707        sage: chebyshev_T(1,z)
708        array([1, 2])
709        sage: chebyshev_T(1,z2)
710        array([[1, 2],
711              [1, 2]])
712        sage: chebyshev_T(1,z3)
713        array([ 1.,  2.,  3.])
714
715        """
716        if isinstance(args[0],numpy.ndarray):
717            raise NotImplementedError("Support for numpy array in \
718            first argument(s) is not supported yet!")
719
720        result = numpy.zeros(args[-1].shape).tolist()
721        if isinstance(args[-1][0],int):
722            for k in xrange(len(args[-1])):
723                result[k] = chebyshev_T(args[0],ZZ(args[-1][k]))
724
725        if isinstance(args[-1][0],float):
726            for k in xrange(len(args[-1])):
727                result[k] = chebyshev_T(args[0],RR(args[-1][k]))
728
729        if isinstance(args[-1][0],numpy.ndarray):
730            for k in xrange(len(args[-1])):
731                result[k] = chebyshev_T(args[0],args[-1][k])
732
733        return numpy.array(result)
734
735    def _derivative_(self, *args, **kwds):
736        """
737        Returns the derivative of chebyshev_T in form of the chebyshev Polynomial
738        of the second kind chebyshev_U
739        EXAMPLES::
740            sage: var('k')
741            k
742            sage: derivative(chebyshev_T(k,x),x)
743            k*chebyshev_U(k - 1, x)
744            sage: derivative(chebyshev_T(3,x),x)
745            12*x^2 - 3
746            sage: derivative(chebyshev_T(k,x),k)
747            Traceback (most recent call last):
748            ...
749            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
750 
751        """
752        diff_param = kwds['diff_param']
753        if diff_param == 0: 
754            raise NotImplementedError(
755"Derivative w.r.t. to the index is not supported, yet, \
756and perhaps never will be...")
757        else:
758            return args[0]*chebyshev_U(args[0]-1,args[1])
759
760   
761chebyshev_T = Func_chebyshev_T()
762
763class Func_chebyshev_U(OrthogonalPolynomial): 
764   
765    """
766    Class for the Chebyshev polynomial of the second kind.
767   
768    REFERENCE:
769
770    - AS, 22.8.3 page 783 and AS 6.1.22 page 256.
771   
772    EXAMPLES::
773   
774        sage: chebyshev_U(3,x)
775        8*x^3 - 4*x
776    """
777
778    def __init__(self):
779        OrthogonalPolynomial.__init__(self,"chebyshev_U",
780nargs = 2,conversions =dict(maxima='chebyshev_u',mathematica='ChebyshevU'))
781       
782    def _clenshaw_method_(self,*args):
783        """
784        Clenshaw method for chebyshev_U (means use the recursion...)
785        This is much faster for numerical evaluation than maxima!
786        See A.S. 227 (p. 782) for details for the recurions
787        """
788        k = args[0]
789        x = args[1]
790           
791        if k == 0:
792            return 1
793        elif k == 1:
794            return 2*x
795        else:
796            help1 = 1
797            help2 = 2*x
798            if is_Expression(x):
799                #raise NotImplementedError("Maxima is faster here!")
800                help3 = 0
801                for j in xrange(0,floor(k/2)+1):
802                    help3 = \
803help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j)/factorial(j)/factorial(k-2*j)
804               
805            else:
806                for j in xrange(0,k-1):
807                    help3 = 2*x*help2 - help1
808                    help1 = help2
809                    help2 = help3
810                   
811            return help3   
812   
813    def _maxima_init_evaled_(self, *args):
814        """
815        Uses 
816        EXAMPLES::
817        sage: x = PolynomialRing(QQ, 'x').gen()
818        sage: chebyshev_T(2,x)
819        2*x^2 - 1
820        """
821        n = args[0]
822        x = args[1]
823        return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
824
825       
826    def _evalf_(self, *args,**kwds):
827        """
828        Evals chebyshev_U
829        numerically with mpmath.
830        EXAMPLES::
831            sage: chebyshev_U(10,3).n(75)
832            4.661117900000000000000e7
833        """
834        try:
835            step_parent = kwds['parent']
836        except KeyError:
837            step_parent = parent(args[-1])
838
839        import sage.libs.mpmath.all as mpmath
840
841        try:
842            precision = step_parent.prec()
843        except AttributeError:
844            precision = mpmath.mp.prec
845
846        return mpmath.call(mpmath.chebyu,args[0],args[-1],prec = precision)
847
848    def _eval_special_values_(self,*args):
849        """
850        Special values known. A.S. 22.4 (p.777).
851        EXAMPLES:
852       
853            sage: var('n')
854            n
855            sage: chebyshev_U(n,1)
856            n + 1
857            sage: chebyshev_U(n,-1)
858            (n + 1)*(-1)^n
859        """
860        if args[-1] == 1:
861            return (args[0]+1)
862       
863        if args[-1] == -1:
864            return (-1)**args[0]*(args[0]+1)
865
866        if (args[-1] == 0):
867            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
868
869        raise ValueError("Value not found")
870
871    def _eval_numpy_(self, *args):
872        #TODO: numpy_eval with help of the new scipy version!!!!
873        #Reason scipy 0.8 supports stable and fast numerical evaluation
874        #of ortho polys
875        #Now this only evaluates the array pointwise, and only the first one...
876        #if isinstance(arg[0], numpy.ndarray).
877        #This is a hack to provide compability!
878        """
879        EXAMPLES::
880        sage: import numpy
881        sage: z = numpy.array([1,2])
882        sage: z2 = numpy.array([[1,2],[1,2]])
883        sage: z3 = numpy.array([1,2,3.])
884        sage: chebyshev_U(1,z)
885        array([2, 4])
886        sage: chebyshev_U(1,z2)
887        array([[2, 4],
888              [2, 4]])
889        sage: chebyshev_U(1,z3)
890        array([ 2.,  4.,  6.])
891
892        """
893        if isinstance(args[0],numpy.ndarray):
894            raise NotImplementedError(
895           "Support for numpy array in first argument(s) is not supported yet!")
896
897        result = numpy.zeros(args[-1].shape).tolist()
898        if isinstance(args[-1][0],int):
899            for k in xrange(len(args[-1])):
900                result[k] = chebyshev_U(args[0],ZZ(args[-1][k]))
901
902        if isinstance(args[-1][0],float):
903            for k in xrange(len(args[-1])):
904                result[k] = chebyshev_U(args[0],RR(args[-1][k]))
905
906        if isinstance(args[-1][0],numpy.ndarray):
907            for k in xrange(len(args[-1])):
908                result[k] = chebyshev_U(args[0],args[-1][k])
909
910        return numpy.array(result)
911
912
913    def _derivative_(self, *args, **kwds):
914        """
915        Returns the derivative of chebyshev_U in form of the chebyshev
916        Polynomials of the first and second kind
917       
918        EXAMPLES::
919            sage: var('k')
920            k
921            sage: derivative(chebyshev_U(k,x),x)
922            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
923            sage: derivative(chebyshev_U(3,x),x)
924            24*x^2 - 4
925            sage: derivative(chebyshev_U(k,x),k)
926            Traceback (most recent call last):
927            ...
928            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
929 
930        """
931        diff_param = kwds['diff_param']
932        if diff_param == 0: 
933            raise NotImplementedError(
934"Derivative w.r.t. to the index is not supported, \
935yet, and perhaps never will be...")
936        else:
937            return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
938                    chebyshev_U(args[0],args[1]))/(args[1]**2-1)
939       
940chebyshev_U = Func_chebyshev_U()
941
942class Func_gegenbauer(OrthogonalPolynomial):
943    """
944    Returns the ultraspherical (or Gegenbauer) polynomial.
945   
946    REFERENCE:
947
948    - AS 22.5.27
949   
950    EXAMPLES::
951   
952        sage: x = PolynomialRing(QQ, 'x').gen()
953        sage: ultraspherical(2,3/2,x)
954        15/2*x^2 - 3/2
955        sage: ultraspherical(2,1/2,x)
956        3/2*x^2 - 1/2
957        sage: ultraspherical(1,1,x)
958        2*x     
959        sage: t = PolynomialRing(RationalField(),"t").gen()
960        sage: gegenbauer(3,2,t)
961        32*t^3 - 12*t
962    """   
963    def __init__(self):
964        OrthogonalPolynomial.__init__(self,"gegenbauer",nargs = 3,
965        conversions =dict(maxima='ultraspherical',mathematica='GegenbauerC'))
966
967    def _clenshaw_method_(self,*args):
968        """
969        Clenshaw method for gegenbauer poly (means use the recursion...)
970        This is much faster for numerical evaluation than maxima!
971        See A.S. 227 (p. 782) for details for the recurions
972        """
973        k = args[0]
974        x = args[-1]
975        alpha = args[1]
976
977        if is_Expression(alpha) or abs(alpha) > numpy.finfo(float).eps:
978            alpha_zero = False
979        else:
980            alpha_zero = True
981
982        if k == 0:
983            return 1
984        elif k == 1:
985            if not alpha_zero:
986                return 2*x*alpha
987            else:
988                return 2*x # It seems that maxima evals this the wrong way!
989                           #(see A.S. 22.4 (p.777))
990        else:
991            help1 = 1
992            if alpha_zero:
993                help2 = 2*x
994                gamma_alpha = 1
995            else:
996                help2 = 2*x*alpha
997                gamma_alpha = gamma(alpha)
998
999            help3 = 0
1000            if is_Expression(x):
1001                for j in xrange(0,floor(k/2)+1):
1002                    help3 = \
1003help3 + (-1)**j*gamma(alpha+k-j)/factorial(j)/factorial(k-2*j)*(2*x)**(k-2*j)/gamma_alpha
1004               
1005            else:
1006                for j in xrange(1,k):
1007                    help3 = 2*(j+alpha)*x*help2 - (j+2*alpha-1)*help1
1008                    help3 = help3/(j+1)
1009                    help1 = help2
1010                    help2 = help3
1011                   
1012            return help3   
1013
1014    def _maxima_init_evaled_(self, *args):
1015        """
1016        Uses 
1017        EXAMPLES::
1018        sage: x = PolynomialRing(QQ, 'x').gen()
1019        sage: ultraspherical(2,1/2,x)
1020        3/2*x^2 - 1/2
1021        """
1022        n = args[0]
1023        a = args[1]
1024        x = args[2]
1025        return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)),\
1026                         locals={'x':x})
1027
1028
1029    def _evalf_(self, *args,**kwds):
1030        """
1031        Evals gegenbauer
1032        numerically with mpmath.
1033        EXAMPLES::
1034            sage: gegenbauer(10,2,3.).n(54)
1035            5.25360702000000e8
1036        """
1037        if args[1] == 0:
1038            raise ValueError("mpmath don't handle alpha = 0 correctly!")
1039           
1040        try:
1041            step_parent = kwds['parent']
1042        except KeyError:
1043            step_parent = parent(args[-1])
1044
1045        import sage.libs.mpmath.all as mpmath
1046
1047        try:
1048            precision = step_parent.prec()
1049        except AttributeError:
1050            precision = mpmath.mp.prec
1051
1052        return mpmath.call(
1053            mpmath.gegenbauer,args[0],args[1],args[-1],prec = precision)
1054
1055    def _eval_special_values_(self,*args):
1056        """
1057        Special values known. A.S. 22.4 (p.777)
1058        EXAMPLES:
1059       
1060            sage: var('n a')
1061            (n, a)
1062            sage: gegenbauer(n,1/2,x)
1063            legendre_P(n, x)
1064            sage: gegenbauer(n,0,x)
1065            1/2*n*chebyshev_T(n, x)
1066            sage: gegenbauer(n,1,x)
1067            chebyshev_U(n, x)
1068            sage: gegenbauer(n,a,1)
1069            binomial(2*a + n - 1, n)
1070            sage: gegenbauer(n,a,-1)
1071            (-1)^n*binomial(2*a + n - 1, n)
1072            sage: gegenbauer(n,a,0)
1073            1/2*((-1)^n + 1)*(-1)^(1/2*n)*gamma(a + 1/2*n)/(gamma(a)*gamma(1/2*n))
1074        """
1075        if args[1] == 0 and args[0] != 0:
1076            return args[0]*chebyshev_T(args[0],args[-1])/2
1077
1078        if args[1] == 0.5:
1079            return legendre_P(args[0],args[-1])
1080
1081        if args[1] == 1:
1082            return chebyshev_U(args[0],args[-1])
1083
1084        if args[-1] == 1:
1085            return binomial(args[0] + 2*args[1] - 1,args[0])
1086       
1087        if args[-1] == -1:
1088            return (-1)**args[0]*binomial(args[0] + 2*args[1] - 1,args[0])
1089
1090        if (args[-1] == 0):
1091            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2*\
1092            gamma(args[1]+args[0]/2)/gamma(args[1])/gamma(args[0]/2)
1093
1094        raise ValueError("Value not found")
1095
1096    def _eval_numpy_(self, *args):
1097        #TODO: numpy_eval with help of the new scipy version!!!!
1098        #Reason scipy 0.8 supports stable and fast numerical evaluation
1099        #of ortho polys
1100        #Now this only evaluates the array pointwise, and only the first one...
1101        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
1102        """
1103        EXAMPLES::
1104        sage: import numpy
1105        sage: z = numpy.array([1,2])
1106        sage: z2 = numpy.array([[1,2],[1,2]])
1107        sage: z3 = numpy.array([1,2,3.])
1108        sage: gegenbauer(2,0,z)
1109        array([1, 7])
1110        sage: gegenbauer(2,1,z)
1111        array([ 3, 15])
1112        sage: gegenbauer(2,1,z2)
1113        array([[ 3, 15],
1114              [ 3, 15]])
1115        sage: gegenbauer(2,1,z3)
1116        array([  3.,  15.,  35.])
1117        """
1118
1119        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
1120            raise NotImplementedError(
1121                "Support for numpy array in first argument(s) is not supported yet!")
1122       
1123        result = numpy.zeros(args[-1].shape).tolist()
1124        if isinstance(args[-1][0],int):
1125            for k in xrange(len(args[-1])):
1126                result[k] = gegenbauer(args[0],args[1],ZZ(args[-1][k]))
1127
1128        if isinstance(args[-1][0],float):
1129            for k in xrange(len(args[-1])):
1130                result[k] = gegenbauer(args[0],args[1],RR(args[-1][k]))
1131
1132        if isinstance(args[-1][0],numpy.ndarray):
1133            for k in xrange(len(args[-1])):
1134                result[k] = gegenbauer(args[0],args[1],args[-1][k])
1135
1136        return numpy.array(result)
1137
1138    def _derivative_(self, *args, **kwds):
1139        """
1140        Returns the derivative of chebyshev_U in form of the chebyshev
1141        Polynomials of the first and second kind
1142       
1143        EXAMPLES::
1144            sage: var('k a')
1145            (k, a)
1146            sage: derivative(gegenbauer(k,a,x),x)
1147            (k*x*gegenbauer(k, a, x) - (2*a + k - 1)*gegenbauer(k - 1, a, x))/(x^2 - 1)
1148            sage: derivative(gegenbauer(4,3,x),x)
1149            960*x^3 - 240*x
1150            sage: derivative(gegenbauer(k,a,x),a)
1151            Traceback (most recent call last):
1152            ...
1153            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
1154 
1155        """
1156        diff_param = kwds['diff_param']
1157        if diff_param in [0,1]: 
1158            raise NotImplementedError(
1159"Derivative w.r.t. to the index is not supported, yet,\
1160and perhaps never will be...")
1161        else:
1162            return (-args[0]*args[-1]*gegenbauer(args[0],args[1],args[2])+\
1163(args[0] + 2*args[1]-1)*gegenbauer(args[0]-1,args[1],args[2]))/(1-args[-1]**2)
1164
1165gegenbauer = Func_gegenbauer()
1166ultraspherical = Func_gegenbauer()
1167
1168class Func_gen_laguerre(OrthogonalPolynomial):
1169    """
1170    Returns the generalized Laguerre polynomial for integers `n > -1`.
1171    Typically, a = 1/2 or a = -1/2.
1172   
1173    REFERENCE:
1174
1175    - table on page 789 in AS.
1176   
1177    EXAMPLES::
1178   
1179        sage: x = PolynomialRing(QQ, 'x').gen()
1180        sage: gen_laguerre(2,1,x)
1181        1/2*x^2 - 3*x + 3
1182        sage: gen_laguerre(2,1/2,x)
1183        1/2*x^2 - 5/2*x + 15/8
1184        sage: gen_laguerre(2,-1/2,x)
1185        1/2*x^2 - 3/2*x + 3/8
1186        sage: gen_laguerre(2,0,x)
1187        1/2*x^2 - 2*x + 1
1188        sage: gen_laguerre(3,0,x)
1189        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1190    """
1191    def __init__(self):
1192        OrthogonalPolynomial.__init__(self,"gen_laguerre",nargs = 3,
1193conversions =dict(maxima='gen_laguerre',mathematica='LaguerreL'))
1194       
1195    def _clenshaw_method_(self,*args):
1196        """
1197        Clenshaw method for gen_laguerre polynomial (means use the recursion...)
1198        This is much faster for numerical evaluation than maxima!
1199        See A.S. 227 (p. 782) for details for the recurions.
1200        Warning: The clanshaw method for the laguerre Polynomials
1201        should only used for exact data types, when high orders are
1202        used, due to weak instabilities of the recursion!
1203        """
1204        k = args[0]
1205        a = args[1]
1206        x = args[-1]
1207       
1208        if k == 0:
1209            return 1
1210        elif k == 1:
1211            return 1-x + a
1212        else:
1213            help1 = 1
1214            help2 = 1-x + a
1215            if is_Expression(x): 
1216                help3 = 0
1217                for j in xrange(0,k+1):
1218                    help3 = help3 + (-1)**j*binomial(k+a,k-j)/factorial(j)*x**j
1219            else:
1220                for j in xrange(1,k):
1221                    help3 = -x*help2 + (2*j+1+a)*help2 - (j+a)*help1
1222                    help3 = help3/(j+1)
1223                    help1 = help2
1224                    help2 = help3
1225               
1226            return help3   
1227
1228    def _maxima_init_evaled_(self, *args):
1229        n = args[0]
1230        a = args[1]
1231        x = args[-1]
1232       
1233        from sage.functions.all import sqrt
1234        _init()
1235        return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)),\
1236                         locals={'x':x})
1237
1238    def _evalf_(self, *args,**kwds):
1239        """
1240        Evals laguerre polynomial
1241        numerically with mpmath.
1242        EXAMPLES::
1243            sage: laguerre(3,5.).n(53)
1244            2.66666666666667
1245        """
1246        try:
1247            step_parent = kwds['parent']
1248        except KeyError:
1249            step_parent = parent(args[-1])
1250
1251        import sage.libs.mpmath.all as mpmath
1252
1253        try:
1254            precision = step_parent.prec()
1255        except AttributeError:
1256            precision = mpmath.mp.prec
1257
1258        return mpmath.call(
1259            mpmath.laguerre,args[0],args[1],args[-1],prec = precision)
1260   
1261    def _eval_special_values_(self,*args):
1262        """
1263        Special values known.
1264        EXAMPLES:
1265           
1266            sage: var('n')
1267            n
1268            sage: gen_laguerre(n,1,0)
1269            binomial(n + 1, n)
1270        """
1271       
1272        if (args[-1] == 0):
1273            try:
1274                return binomial(args[0]+args[1],args[0])
1275            except TypeError:
1276                pass
1277
1278        raise ValueError("Value not found")
1279
1280    def _eval_numpy_(self, *args):
1281        #TODO: numpy_eval with help of the new scipy version!!!!
1282        #Reason scipy 0.8 supports stable and fast numerical evaluation
1283        #of ortho polys
1284        #Now this only evaluates the array pointwise, and only the first one...
1285        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
1286        """
1287        EXAMPLES::
1288        sage: import numpy
1289        sage: z = numpy.array([1,2])
1290        sage: z2 = numpy.array([[1,2],[1,2]])
1291        sage: z3 = numpy.array([1,2,3.])
1292        sage: gen_laguerre(1,1,z)
1293        array([1, 0])
1294        sage: gen_laguerre(1,1,z2)
1295        array([[1, 0],
1296              [1, 0]])
1297        sage: gen_laguerre(1,1,z3)
1298        array([ 1.,  0., -1.])
1299        """
1300
1301        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
1302            raise NotImplementedError(
1303           "Support for numpy array in first argument(s) is not supported yet!")
1304       
1305        result = numpy.zeros(args[-1].shape).tolist()
1306        if isinstance(args[-1][0],int):
1307            for k in xrange(len(args[-1])):
1308                result[k] = gen_laguerre(args[0],args[1],ZZ(args[-1][k]))
1309
1310        if isinstance(args[-1][0],float):
1311            for k in xrange(len(args[-1])):
1312                result[k] = gen_laguerre(args[0],args[1],RR(args[-1][k]))
1313
1314        if isinstance(args[-1][0],numpy.ndarray):
1315            for k in xrange(len(args[-1])):
1316                result[k] = gen_laguerre(args[0],args[1],args[-1][k])
1317
1318        return numpy.array(result)
1319
1320   
1321
1322    def _derivative_(self,*args,**kwds):
1323        """return the derivative of gen_laguerre in
1324        form of the gen. Laguerre Polynomial.
1325        EXAMPLES::
1326        sage: n = var('n')
1327        sage: derivative(gen_laguerre(3,1,x),x)
1328        -1/2*x^2 + 4*x - 6
1329        sage: derivative(gen_laguerre(n,1,x),x)
1330        -((n + 1)*gen_laguerre(n - 1, 1, x) - n*gen_laguerre(n, 1, x))/x
1331        """
1332        diff_param = kwds['diff_param']
1333        if diff_param == 0: 
1334            raise NotImplementedError(
1335"Derivative w.r.t. to the index is not supported, \
1336yet, and perhaps never will be...")
1337        else:
1338            return (args[0]*gen_laguerre(args[0],args[1],args[-1])-(args[0]+args[1])*\
1339                    gen_laguerre(args[0]-1,args[1],args[-1]))/args[-1]
1340
1341gen_laguerre = Func_gen_laguerre()
1342
1343class Func_hermite(OrthogonalPolynomial):
1344    """
1345    Returns the Hermite polynomial for integers `n > -1`.
1346   
1347    REFERENCE:
1348
1349    - AS 22.5.40 and 22.5.41, page 779.
1350   
1351    #sage: S.<y> = PolynomialRing(RR)
1352    #   <- Here a bug of the base class BuiltinFunction occours!
1353
1354    EXAMPLES::
1355   
1356        sage: x = PolynomialRing(QQ, 'x').gen()
1357        sage: hermite(2,x)
1358        4*x^2 - 2
1359        sage: hermite(3,x)
1360        8*x^3 - 12*x
1361        sage: hermite(3,2)
1362        40
1363        sage: y = var('y')
1364        sage: hermite(3,y).polynomial(PolynomialRing(RR, 'y'))
1365        8.00000000000000*y^3 - 12.0000000000000*y
1366        sage: R.<x,y> = QQ[]
1367        sage: hermite(3,y^2)
1368        8*y^6 - 12*y^2
1369        sage: w = var('w')
1370        sage: hermite(3,2*w)
1371        64*w^3 - 24*w
1372    """
1373   
1374    def __init__(self):
1375        OrthogonalPolynomial.__init__(self,"hermite",nargs = 2,
1376conversions =dict(maxima='hermite',mathematica='HermiteH'))
1377
1378    def _clenshaw_method_(self,*args):
1379        """
1380        Clenshaw method for hermite polynomial (means use the recursion...)
1381        See A.S. 227 (p. 782) for details for the recurions
1382        For the symbolic evaluation, maxima seems to be quite fast.
1383        The break even point between the recursion and Maxima is about
1384        n = 25
1385        """
1386        k = args[0]
1387        x = args[1]
1388           
1389        if k == 0:
1390            return 1
1391        elif k == 1:
1392            return 2*x
1393        else:
1394            help1 = 1
1395            help2 = 2*x
1396            if is_Expression(x):
1397                help3 = 0
1398                for j in xrange(0,floor(k/2)+1):
1399                    help3 = help3 + (-1)**j*(2*x)**(k-2*j)/factorial(j)/\
1400                            factorial(k-2*j)
1401                help3 = help3*factorial(k)
1402            else:
1403                for j in xrange(1,k):
1404                    help3 = 2*x*help2 - 2*j*help1
1405                    help1 = help2
1406                    help2 = help3
1407                   
1408            return help3   
1409   
1410           
1411    def _evalf_(self, *args,**kwds):
1412        """
1413        Evals hermite
1414        numerically with mpmath.
1415        EXAMPLES::
1416            sage: hermite(3,2.).n(74)
1417            40.0000000000000000000
1418        """
1419        try:
1420            step_parent = kwds['parent']
1421        except KeyError:
1422            step_parent = parent(args[-1])
1423
1424        import sage.libs.mpmath.all as mpmath
1425
1426        try:
1427            precision = step_parent.prec()
1428        except AttributeError:
1429            precision = mpmath.mp.prec
1430
1431        return mpmath.call(mpmath.hermite,args[0],args[-1],prec = precision)
1432   
1433    def _eval_special_values_(self,*args):
1434        """
1435        Special values known. A.S. 22.4 (p.777)
1436        EXAMPLES:
1437       
1438            sage: var('n')
1439            n
1440            sage: hermite(n,0)
1441            ((-1)^n + 1)*(-1)^(1/2*n)*factorial(n)/gamma(1/2*n + 1)
1442        """
1443
1444        if (args[-1] == 0):
1445            return (1+(-1)**args[0])*(-1)**(args[0]/2)*\
1446                   factorial(args[0])/gamma(args[0]/2+1)
1447
1448        raise ValueError("Value not found")
1449
1450    def _maxima_init_evaled_(self, *args):
1451        """
1452        Old maxima method.
1453        """
1454        n = args[0]
1455        x = args[1]
1456        return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
1457
1458    def _eval_numpy_(self, *args):
1459        #TODO: numpy_eval with help of the new scipy version!!!!
1460        #Reason scipy 0.8 supports stable and fast numerical evaluation
1461        #of ortho polys
1462        #Now this only evaluates the array pointwise, and only the first one...
1463        #if isinstance(arg[0], numpy.ndarray).
1464        #This is a hack to provide compability!
1465        """
1466        EXAMPLES::
1467        sage: import numpy
1468        sage: z = numpy.array([1,2])
1469        sage: z2 = numpy.array([[1,2],[1,2]])
1470        sage: z3 = numpy.array([1,2,3.])
1471        sage: hermite(1,z)
1472        array([2, 4])
1473        sage: hermite(1,z2)
1474        array([[2, 4],
1475              [2, 4]])
1476        sage: hermite(1,z3)
1477        array([ 2.,  4.,  6.])
1478
1479        """
1480        if isinstance(args[0],numpy.ndarray):
1481            raise NotImplementedError(
1482          "Support for numpy array in first argument(s) is not supported yet!")
1483
1484        result = numpy.zeros(args[-1].shape).tolist()
1485        if isinstance(args[-1][0],int):
1486            for k in xrange(len(args[-1])):
1487                result[k] = hermite(args[0],ZZ(args[-1][k]))
1488
1489        if isinstance(args[-1][0],float):
1490            for k in xrange(len(args[-1])):
1491                result[k] = hermite(args[0],RR(args[-1][k]))
1492
1493        if isinstance(args[-1][0],numpy.ndarray):
1494            for k in xrange(len(args[-1])):
1495                result[k] = hermite(args[0],args[-1][k])
1496
1497        return numpy.array(result)
1498
1499
1500    def _derivative_(self, *args, **kwds):
1501        """
1502        Returns the derivative of the hermite polynomial in form of the chebyshev
1503        Polynomials of the first and second kind
1504       
1505        EXAMPLES::
1506            sage: var('k')
1507            k
1508            sage: derivative(hermite(k,x),x)
1509            2*k*hermite(k - 1, x)
1510             
1511        """
1512        diff_param = kwds['diff_param']
1513        if diff_param == 0: 
1514            raise NotImplementedError(
1515"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1516        else:
1517            return 2*args[0]*hermite(args[0]-1,args[1])
1518 
1519hermite = Func_hermite()       
1520
1521
1522class Func_jacobi_P(OrthogonalPolynomial):
1523    r"""
1524    Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for
1525    integers `n > -1` and a and b symbolic or `a > -1`
1526    and `b > -1`. The Jacobi polynomials are actually defined
1527    for all a and b. However, the Jacobi polynomial weight
1528    `(1-x)^a(1+x)^b` isn't integrable for `a \leq -1`
1529    or `b \leq -1`.
1530   
1531    REFERENCE:
1532
1533    - table on page 789 in AS.
1534   
1535    EXAMPLES::
1536   
1537        sage: x = PolynomialRing(QQ, 'x').gen()
1538        sage: jacobi_P(2,0,0,x)
1539        3/2*x^2 - 1/2
1540        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
1541        5.009999999999998
1542    """
1543
1544    def __init__(self):
1545        OrthogonalPolynomial.__init__(self,"jacobi_P",nargs = 4,
1546conversions =dict(maxima='jacobi_p',mathematica='JacobiP'))
1547   
1548    def _clenshaw_method_(self,*args):
1549        """
1550        Clenshaw method for jacobi_P (means use the recursion,
1551        or sum formula)
1552        This is much faster for numerical evaluation than maxima!
1553        See A.S. 227 (p. 782) for details for the recurions.
1554        Warning: The clanshaw method for the Jacobi Polynomials
1555        should only used for exact data types, when high orders are
1556        used, due to weak instabilities of the recursion!
1557        """ 
1558        k = args[0]
1559        x = args[-1]
1560        alpha = args[1]
1561        beta = args[2]
1562
1563        if k == 0:
1564            return 1
1565        elif k == 1:
1566            return (alpha-beta + (alpha+beta+2)*x)/2
1567        else:
1568           
1569            if is_Expression(x) or is_Expression(alpha) or is_Expression(beta):
1570                #Here we use the sum formula of jacobi_P it seems this is rather
1571                #optimal for use.
1572                help1 = gamma(alpha+k+1)/factorial(k)/gamma(alpha+beta+k+1)
1573                help2 = 0
1574                for j in xrange(0,k+1):
1575                    help2 = help2 + binomial(k,j)*gamma(alpha+beta+k+j+1)/\
1576                            gamma(alpha+j+1)*((x-1)/2)**j
1577                return help1*help2
1578            else:
1579                help1 = 1
1580                help2 = (alpha-beta + (alpha+beta+2)*x)/2
1581
1582                for j in xrange(1,k):
1583                    a1 = 2*(j+1)*(j+alpha+beta+1)*(2*j+alpha+beta)
1584                    a2 = (2*j+alpha+beta+1)*(alpha**2-beta**2)
1585                    a3 = (2*j+alpha+beta)*(2*j+alpha+beta+1)*(2*j+alpha+beta+2)
1586                    a4 = 2*(j+alpha)*(j+beta)*(2*j+alpha+beta+2)
1587                    help3 = (a2+a3*x)*help2 - a4*help1
1588                    help3 = help3/a1
1589                    help1 = help2
1590                    help2 = help3
1591               
1592            return help3   
1593
1594    def _maxima_init_evaled_(self, *args):
1595        """
1596        Old maxima method.
1597        """
1598        n = args[0]
1599        a = args[1]
1600        b = args[2]
1601        x = args[3]
1602        return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)),\
1603                         locals={'x':x})
1604
1605    def _evalf_(self, *args,**kwds):
1606        """
1607        Evals jacobi_P
1608        numerically with mpmath.
1609        EXAMPLES::
1610            sage: jacobi_P(10,2,3,3).n(75)
1611            1.322776620000000000000e8
1612        """
1613        try:
1614            step_parent = kwds['parent']
1615        except KeyError:
1616            step_parent = parent(args[-1])
1617
1618        import sage.libs.mpmath.all as mpmath
1619
1620        try:
1621            precision = step_parent.prec()
1622        except AttributeError:
1623            precision = mpmath.mp.prec
1624
1625        return mpmath.call(mpmath.jacobi,args[0],args[1],args[2],args[-1],\
1626                           prec = precision)
1627
1628    def _eval_special_values_(self,*args):
1629        """
1630        Special values known. A.S. 22.4 (p.777)
1631        EXAMPLES:
1632           
1633            sage: var('n')
1634            n
1635            sage: legendre_P(n,1)
1636            1
1637            sage: legendre_P(n,-1)
1638            (-1)^n
1639        """
1640        if args[-1] == 1:
1641            return binomial(args[0]+args[1],args[0])
1642       
1643        if args[-1] == -1:
1644            return (-1)**args[0]*binomial(args[0]+args[1],args[0])
1645       
1646        if args[1] == 0 and args[2] == 0:
1647            return legendre_P(args[0],args[-1])
1648
1649        if args[1] == -0.5 and args[2] == -0.5:
1650            try:
1651                return binomial(2*args[0],args[0])*\
1652                       chebyshev_T(args[0],args[-1])/4**args[0]
1653            except TypeError:
1654                pass
1655
1656        raise ValueError("Value not found")
1657
1658    def _eval_numpy_(self, *args):
1659        #TODO: numpy_eval with help of the new scipy version!!!!
1660        #Reason scipy 0.8 supports stable and fast numerical evaluation
1661        #of ortho polys
1662        #Now this only evaluates the array pointwise, and only the first one...
1663        #if isinstance(arg[0], numpy.ndarray).
1664        #This is a hack to provide compability!
1665        """
1666        EXAMPLES::
1667        sage: import numpy
1668        sage: z = numpy.array([1,2])
1669        sage: z2 = numpy.array([[1,2],[1,2]])
1670        sage: z3 = numpy.array([1,2,3.])
1671        sage: gen_laguerre(1,1,z)
1672        array([1, 0])
1673        sage: gen_laguerre(1,1,z2)
1674        array([[1, 0],
1675              [1, 0]])
1676        sage: gen_laguerre(1,1,z3)
1677        array([ 1.,  0., -1.])
1678        """
1679
1680        if isinstance(args[0],numpy.ndarray) or \
1681               isinstance(args[1],numpy.ndarray) or \
1682               isinstance(args[2],numpy.ndarray):
1683            raise NotImplementedError(
1684           "Support for numpy array in first argument(s) is not supported yet!")
1685       
1686        result = numpy.zeros(args[-1].shape).tolist()
1687        if isinstance(args[-1][0],int):
1688            for k in xrange(len(args[-1])):
1689                result[k] = jacobi_P(args[0],args[1],args[2],ZZ(args[-1][k]))
1690
1691        if isinstance(args[-1][0],float):
1692            for k in xrange(len(args[-1])):
1693                result[k] = jacobi_P(args[0],args[1],args[2],RR(args[-1][k]))
1694
1695        if isinstance(args[-1][0],numpy.ndarray):
1696            for k in xrange(len(args[-1])):
1697                result[k] = jacobi_P(args[0],args[1],args[2],args[-1][k])
1698
1699        return numpy.array(result)
1700
1701    def _derivative_(self, *args, **kwds):
1702        """
1703        Returns the derivative of jacobi_P in form of jacobi_polynomials
1704       
1705        EXAMPLES::
1706            sage: var('k a b')
1707            (k, a, b)
1708            sage: derivative(jacobi_P(k,a,b,x),x)
1709            -(2*(b + k)*(a + k)*jacobi_P(k - 1, a, b, x) - ((a + b + 2*k)*x - a + b)*k*jacobi_P(k, a, b, x))/((x^2 - 1)*(a + b + 2*k))
1710            sage: derivative(jacobi_P(2,1,3,x),x)
1711            14*x - 7/2
1712            sage: derivative(jacobi_P(k,a,b,x),a)
1713            Traceback (most recent call last):
1714            ...
1715            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
1716 
1717        """
1718        diff_param = kwds['diff_param']
1719        if diff_param in [0,1,2]: 
1720            raise NotImplementedError(
1721"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1722        else:
1723            return (args[0]*(args[1]-args[2]-(2*args[0]+args[1]+args[2])*args[-1])*\
1724                    jacobi_P(args[0],args[1],args[2],args[-1])+ 2*(args[0]+args[1])*\
1725                    (args[0]+args[2])*jacobi_P(args[0]-1,args[1],args[2],args[-1]))/\
1726                    (2*args[0]+args[1]+args[2])/(1-args[-1]**2)
1727       
1728       
1729jacobi_P = Func_jacobi_P()       
1730
1731class Func_laguerre(OrthogonalPolynomial):
1732    """
1733    Returns the Laguerre polynomial for integers `n > -1`.
1734   
1735    REFERENCE:
1736
1737    - AS 22.5.16, page 778 and AS page 789.
1738   
1739    EXAMPLES::
1740   
1741        sage: x = PolynomialRing(QQ, 'x').gen()
1742        sage: laguerre(2,x)
1743        1/2*x^2 - 2*x + 1
1744        sage: laguerre(3,x)
1745        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1746        sage: laguerre(2,2)
1747        -1
1748    """
1749    def __init__(self):
1750        OrthogonalPolynomial.__init__(self,"laguerre",nargs = 2,
1751conversions =dict(maxima='laguerre',mathematica='LaguerreL'))
1752       
1753    def _clenshaw_method_(self,*args):
1754        """
1755        Clenshaw method for laguerre polynomial (means use the recursion...)
1756        This is much faster for numerical evaluation than maxima!
1757        See A.S. 227 (p. 782) for details for the recurions.
1758        Warning: The clanshaw method for the laguerre Polynomials
1759        should only used for exact data types, when high orders are
1760        used, due to weak instabilities of the recursion!
1761        """
1762        k = args[0]
1763        x = args[-1]
1764       
1765        if k == 0:
1766            return 1
1767        elif k == 1:
1768            return 1-x
1769        else:
1770            help1 = 1
1771            help2 = 1-x
1772            if is_Expression(x):
1773                help3 = 0
1774                for j in xrange(0,k+1):
1775                    help3 = help3 + (-1)**j*binomial(k,k-j)/factorial(j)*x**j
1776            else:
1777                for j in xrange(1,k):
1778                    help3 = -x*help2 + (2*j+1)*help2 - j*help1
1779                    help3 = help3/(j+1)
1780                    help1 = help2
1781                    help2 = help3
1782               
1783            return help3   
1784
1785    def _maxima_init_evaled_(self, *args):
1786        n = args[0]
1787        x = args[1]
1788        return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
1789
1790    def _evalf_(self, *args,**kwds):
1791        """
1792        Evals laguerre polynomial
1793        numerically with mpmath.
1794        EXAMPLES::
1795            sage: laguerre(3,5.).n(53)
1796            2.66666666666667
1797        """
1798        try:
1799            step_parent = kwds['parent']
1800        except KeyError:
1801            step_parent = parent(args[-1])
1802
1803        import sage.libs.mpmath.all as mpmath
1804
1805        try:
1806            precision = step_parent.prec()
1807        except AttributeError:
1808            precision = mpmath.mp.prec
1809
1810        return mpmath.call(mpmath.laguerre,args[0],0,args[-1],prec = precision)
1811   
1812    def _eval_special_values_(self,*args):
1813        """
1814        Special values known.
1815        EXAMPLES:
1816           
1817            sage: var('n')
1818            n
1819            sage: laguerre(n,0)
1820            1
1821        """
1822       
1823        if (args[-1] == 0):
1824            try:
1825                return 1
1826            except TypeError:
1827                pass
1828
1829        raise ValueError("Value not found")
1830 
1831
1832    def _eval_numpy_(self, *args):
1833        #TODO: numpy_eval with help of the new scipy version!!!!
1834        #Reason scipy 0.8 supports stable and fast numerical evaluation
1835        #of ortho polys
1836        #Now this only evaluates the array pointwise, and only the first one...
1837        #if isinstance(arg[0], numpy.ndarray).
1838        #This is a hack to provide compability!
1839        """
1840        EXAMPLES::
1841        sage: import numpy
1842        sage: z = numpy.array([1,2])
1843        sage: z2 = numpy.array([[1,2],[1,2]])
1844        sage: z3 = numpy.array([1,2,3.])
1845        sage: laguerre(1,z)
1846        array([ 0, -1])
1847        sage: laguerre(1,z2)
1848        array([[ 0, -1],
1849               [ 0, -1]])
1850        sage: laguerre(1,z3)
1851        array([ 0., -1., -2.])
1852
1853        """
1854        if isinstance(args[0],numpy.ndarray):
1855            raise NotImplementedError(
1856           "Support for numpy array in first argument(s) is not supported yet!")
1857
1858        result = numpy.zeros(args[-1].shape).tolist()
1859        if isinstance(args[-1][0],int):
1860            for k in xrange(len(args[-1])):
1861                result[k] = laguerre(args[0],ZZ(args[-1][k]))
1862
1863        if isinstance(args[-1][0],float):
1864            for k in xrange(len(args[-1])):
1865                result[k] = laguerre(args[0],RR(args[-1][k]))
1866
1867        if isinstance(args[-1][0],numpy.ndarray):
1868            for k in xrange(len(args[-1])):
1869                result[k] = laguerre(args[0],args[-1][k])
1870
1871        return numpy.array(result)
1872
1873    def _derivative_(self,*args,**kwds):
1874        """return the derivative of laguerre in
1875        form of the Laguerre Polynomial.
1876        EXAMPLES::
1877        sage: n = var('n')
1878        sage: derivative(laguerre(3,x),x)
1879        -1/2*x^2 + 3*x - 3
1880        sage: derivative(laguerre(n,x),x)
1881        -(n*laguerre(n - 1, x) - n*laguerre(n, x))/x
1882        """
1883        diff_param = kwds['diff_param']
1884        if diff_param == 0: 
1885            raise NotImplementedError(
1886            "Derivative w.r.t. to the index is not supported, \
1887             yet, and perhaps never will be...")
1888        else:
1889            return (args[0]*laguerre(args[0],args[-1])-args[0]*\
1890                    laguerre(args[0]-1,args[-1]))/args[1]
1891
1892laguerre = Func_laguerre()
1893
1894class Func_legendre_P(OrthogonalPolynomial):
1895    """
1896    Returns the Legendre polynomial of the first kind..
1897   
1898    REFERENCE:
1899
1900    - AS 22.5.35 page 779.
1901   
1902    EXAMPLES::
1903   
1904        sage: P.<t> = QQ[]
1905        sage: legendre_P(2,t)
1906        3/2*t^2 - 1/2
1907        sage: legendre_P(3, 1.1)
1908        1.67750000000000
1909        sage: legendre_P(3, 1 + I)
1910        7/2*I - 13/2
1911        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
1912        Traceback (most recent call last):
1913        ...
1914        TypeError: cannot coerce arguments: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Integer Ring to Symbolic Ring
1915        sage: legendre_P(3, GF(11)(5))
1916        8
1917    """
1918    def __init__(self):
1919        OrthogonalPolynomial.__init__(self,"legendre_P",nargs = 2,
1920conversions =dict(maxima='legendre_p',mathematica='LegendreP'))
1921
1922    def _clenshaw_method_(self,*args):
1923        """
1924        Clenshaw method for legendre_P (means use the recursion...)
1925        This is much faster for numerical evaluation than maxima!
1926        See A.S. 227 (p. 782) for details for the recurions.
1927        Warning: The clanshaw method for the Legendre Polynomials
1928        should only used for exact data types, when high orders are
1929        used, due to weak instabilities of the recursion!
1930        """
1931        k = args[0]
1932        x = args[-1]
1933       
1934        if k == 0:
1935            return 1
1936        elif k == 1:
1937            return x
1938        else:
1939            help1 = 1
1940            help2 = x
1941            if is_Expression(x):
1942                #raise NotImplementedError("Maxima is faster here...")
1943                help1 = ZZ(2**k) #Workarround because of segmentation fault...
1944                help3 = 0
1945                for j in xrange(0,floor(k/2)+1):
1946                    help3 = help3 + (-1)**j*x**(k-2*j)*binomial(k,j)*\
1947                            binomial(2*(k-j),k)
1948               
1949                help3 = help3/help1
1950            else:
1951                for j in xrange(1,k):
1952                    help3 = (2*j+1)*x*help2 - j*help1
1953                    help3 = help3/(j+1)
1954                    help1 = help2
1955                    help2 = help3
1956               
1957            return help3   
1958
1959    def _maxima_init_evaled_(self, *args):
1960        n = args[0]
1961        x = args[1]
1962        return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)),\
1963                         locals={'x':x})
1964
1965    def _evalf_(self, *args,**kwds):
1966        """
1967        Evals legendre_P
1968        numerically with mpmath.
1969        EXAMPLES::
1970            sage: legendre_P(10,3).n(75)
1971            8.097453000000000000000e6
1972        """
1973        try:
1974            step_parent = kwds['parent']
1975        except KeyError:
1976            step_parent = parent(args[-1])
1977
1978        import sage.libs.mpmath.all as mpmath
1979
1980        try:
1981            precision = step_parent.prec()
1982        except AttributeError:
1983            precision = mpmath.mp.prec
1984
1985        return mpmath.call(
1986            mpmath.legendre,args[0],args[-1],prec = precision)
1987   
1988
1989    def _eval_special_values_(self,*args):
1990        """
1991        Special values known.
1992        EXAMPLES:
1993           
1994            sage: var('n')
1995            n
1996            sage: legendre_P(n,1)
1997            1
1998            sage: legendre_P(n,-1)
1999            (-1)^n
2000        """
2001        if args[-1] == 1:
2002            return 1
2003       
2004        if args[-1] == -1:
2005            return (-1)**args[0]
2006       
2007        if (args[-1] == 0):
2008            try:
2009                return (1+(-1)**args[0])/2*binomial(args[0],args[0]/2)/4**(args[0]/2)
2010            except TypeError:
2011                pass
2012
2013        raise ValueError("Value not found")
2014
2015    def _eval_numpy_(self, *args):
2016        #TODO: numpy_eval with help of the new scipy version!!!!
2017        #Reason scipy 0.8 supports stable and fast numerical evaluation
2018        #of ortho polys
2019        #Now this only evaluates the array pointwise, and only the first one...
2020        #if isinstance(arg[0], numpy.ndarray).
2021        #This is a hack to provide compability!
2022        """
2023        EXAMPLES::
2024        sage: import numpy
2025        sage: z = numpy.array([1,2])
2026        sage: z2 = numpy.array([[1,2],[1,2]])
2027        sage: z3 = numpy.array([1,2,3.])
2028        sage: legendre_P(1,z)
2029        array([1, 2])
2030        sage: legendre_P(1,z2)
2031        array([[1, 2],
2032              [1, 2]])
2033        sage: legendre_P(1,z3)
2034        array([ 1.,  2.,  3.])
2035
2036        """
2037        if isinstance(args[0],numpy.ndarray):
2038            raise NotImplementedError(
2039           "Support for numpy array in first argument(s) is not supported yet!")
2040
2041        result = numpy.zeros(args[-1].shape).tolist()
2042        if isinstance(args[-1][0],int):
2043            for k in xrange(len(args[-1])):
2044                result[k] = legendre_P(args[0],ZZ(args[-1][k]))
2045
2046        if isinstance(args[-1][0],float):
2047            for k in xrange(len(args[-1])):
2048                result[k] = legendre_P(args[0],RR(args[-1][k]))
2049
2050        if isinstance(args[-1][0],numpy.ndarray):
2051            for k in xrange(len(args[-1])):
2052                result[k] = legendre_P(args[0],args[-1][k])
2053
2054        return numpy.array(result)
2055
2056    def _derivative_(self,*args,**kwds):
2057        """return the derivative of legendre_P in
2058        form of the Legendre Polynomial.
2059        EXAMPLES::
2060        sage: n = var('n')
2061        sage: derivative(legendre_P(n,x),x)
2062        (n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x^2 - 1)
2063        sage: derivative(legendre_P(3,x),x)
2064        15/2*x^2 - 3/2
2065        """
2066        diff_param = kwds['diff_param']
2067        if diff_param == 0: 
2068            raise NotImplementedError(
2069"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
2070        else:
2071            return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*\
2072                    legendre_P(args[0],args[-1]))/(1-args[-1]**2)
2073       
2074
2075legendre_P = Func_legendre_P()     
2076
2077class Func_legendre_Q(OrthogonalPolynomial):
2078    """
2079    Return the chebyshev function of the second kind.
2080   
2081    REFERENCE:
2082   
2083    - A.S. 8 (p. 332)
2084   
2085    EXAMPLES::
2086   
2087        sage: P.<t> = QQ[]
2088        sage: legendre_Q(2, t)
2089        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
2090        sage: legendre_Q(3, 0.5)
2091        -0.198654771479482
2092        sage: legendre_Q(4, 2)
2093        443/16*I*pi + 443/16*log(3) - 365/12
2094        sage: legendre_Q(4, 2.0)
2095        0.00116107583162041 + 86.9828465962674*I
2096
2097    """
2098    def __init__(self):
2099        OrthogonalPolynomial.__init__(self,"legendre_Q",nargs = 2,
2100conversions =dict(maxima='legendre_q',mathematica='LegendreQ'))
2101
2102    def _maxima_init_evaled_(self, *args):
2103        """
2104        Maxima seems just fine for legendre Q. So we use it here!
2105        """
2106        n = args[0]
2107        x = args[1]
2108        return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)),\
2109                         locals={'x':x})
2110
2111    def _clenshaw_method_(self,*args):
2112        """
2113        Clenshaw method for legendre_q (means use the recursion...)
2114        This is much faster for numerical evaluation than maxima!
2115        See A.S. 8.5.3 (p. 334) for details for the recurions.
2116        Warning: The clanshaw method for the Legendre fUNCTIONS
2117        should only used for exact data types, when high orders are
2118        used, due to weak instabilities of the recursion!
2119        """
2120        raise NotImplementedError("Function not ready yet...")
2121       
2122        k = args[0]
2123        x = args[-1]
2124       
2125        if k == 0:
2126            return ln((1+x)/(1-x))/2
2127        elif k == 1:
2128            return x/2*ln((1+x)/(1-x))-1
2129        else:
2130            if is_Expression(x):
2131                raise NotImplementedError("Maxima works fine here!") 
2132                #it seems that the old method just works fine here...
2133                #raise NotImplementedError("clenshaw does not work well...")
2134            else:
2135                help1 = ln((1+x)/(1-x))/2
2136                help2 = x/2*ln((1+x)/(1-x))-1
2137
2138                for j in xrange(1,k):
2139                    help3 = (2*j+1)*x*help2 - j*help1
2140                    help3 = help3/(j+1)
2141                    help1 = help2
2142                    help2 = help3
2143               
2144            return help3   
2145
2146    def _eval_special_values_(self,*args):
2147        """
2148        Special values known.
2149        EXAMPLES:
2150           
2151            sage: var('n')
2152            n
2153            sage: legendre_Q(n,0)
2154            -1/2*sqrt(pi)*gamma(1/2*n + 1/2)*sin(1/2*pi*n)/gamma(1/2*n + 1)
2155        """
2156        if args[-1] == 1:
2157            return NaN
2158       
2159        if args[-1] == -1:
2160            return NaN
2161       
2162        if (args[-1] == 0):
2163            if is_Expression(args[0]): 
2164                try:
2165                    return -(sqrt(SR.pi()))/2*sin(SR.pi()/2*args[0])*\
2166                           gamma((args[0]+1)/2)/gamma(args[0]/2 + 1)
2167                except TypeError:
2168                    pass
2169            else:
2170                return -(sqrt(math.pi))/2*sin(math.pi/2*args[0])*\
2171                       gamma((args[0]+1)/2)/gamma(args[0]/2. + 1)
2172
2173        raise ValueError("Value not found")
2174
2175    def _evalf_(self, *args,**kwds):
2176        """
2177        Evals legendre_Q
2178        numerically with mpmath.
2179        EXAMPLES::
2180            sage: legendre_Q(10,3).n(75)
2181            2.079454941572578263731e-9 + 1.271944942879431601408e7*I
2182        """
2183
2184        try:
2185            step_parent = kwds['parent']
2186        except KeyError:
2187            step_parent = parent(args[-1])
2188
2189        import sage.libs.mpmath.all as mpmath
2190
2191        try:
2192            precision = step_parent.prec()
2193        except AttributeError:
2194            precision = mpmath.mp.prec
2195
2196        return conjugate(
2197            mpmath.call(mpmath.legenq,args[0],0,args[-1],prec = precision))
2198        #it seems that mpmath uses here a different branch of the logarithm
2199
2200    def _eval_numpy_(self, *args):
2201        #TODO: numpy_eval with help of the new scipy version!!!!
2202        #Reason scipy 0.8 supports stable and fast numerical evaluation
2203        #of ortho polys
2204        #Now this only evaluates the array pointwise, and only the first one...
2205        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
2206        """
2207        EXAMPLES::
2208        sage: import numpy
2209        sage: z = numpy.array([1,2])
2210        sage: z2 = numpy.array([[1,2],[1,2]])
2211        sage: z3 = numpy.array([1,2,3.])
2212        sage: legendre_Q(1,z/5.)
2213        array([-0.95945349, -0.83054043])
2214        sage: legendre_Q(1,z2/5.)
2215        array([[-0.95945349, -0.83054043],
2216               [-0.95945349, -0.83054043]])
2217        sage: legendre_Q(1,z3/5.)
2218        array([-0.95945349, -0.83054043, -0.58411169])
2219
2220        """
2221        if isinstance(args[0],numpy.ndarray):
2222            raise NotImplementedError(
2223           "Support for numpy array in first argument(s) is not supported yet!")
2224
2225        result = numpy.zeros(args[-1].shape).tolist()
2226        if isinstance(args[-1][0],int):
2227            for k in xrange(len(args[-1])):
2228                result[k] = legendre_Q(args[0],ZZ(args[-1][k]))
2229
2230        if isinstance(args[-1][0],float):
2231            for k in xrange(len(args[-1])):
2232                result[k] = legendre_Q(args[0],RR(args[-1][k]))
2233
2234        if isinstance(args[-1][0],numpy.ndarray):
2235            for k in xrange(len(args[-1])):
2236                result[k] = legendre_Q(args[0],args[-1][k])
2237
2238        return numpy.array(result)
2239
2240    def _derivative_(self,*args,**kwds):
2241        """return the derivative of legendre_Q in
2242        form of the Legendre Function.
2243        EXAMPLES::
2244        n = var('n')
2245        derivative(legendre_Q(n,x),x)
2246        (n*x*legendre_Q(n, x) - n*legendre_Q(n - 1, x))/(x^2 - 1)
2247        sage: derivative(legendre_Q(3,x),x)
2248        5/4*(x - 1)*(1/(x - 1) - (x + 1)/(x - 1)^2)*x^3/(x + 1) + 15/4*x^2*log(-(x + 1)/(x - 1)) - 3/4*(x - 1)*(1/(x - 1) - (x + 1)/(x - 1)^2)*x/(x + 1) - 5*x - 3/4*log(-(x + 1)/(x - 1))
2249        """
2250        diff_param = kwds['diff_param']
2251        if diff_param == 0: 
2252            raise NotImplementedError(
2253"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
2254        else:
2255            return (args[0]*args[-1]*legendre_Q(args[0],args[-1])-args[0]*\
2256                    legendre_Q(args[0]-1,args[-1]))/(args[-1]**2-1)
2257
2258legendre_Q = Func_legendre_Q()
2259
2260
2261class Func_gen_legendre_P(OrthogonalPolynomial):
2262
2263    def __init__(self):
2264        OrthogonalPolynomial.__init__(self,"gen_legendre_P",nargs = 3,
2265conversions =dict(maxima='assoc_legendre_p',mathematica='LegendreP'))
2266
2267    def _evalf_(self, *args,**kwds):
2268        """
2269        Evals gen_legendre_P
2270        numerically with mpmath.
2271        EXAMPLES::
2272            sage: gen_legendre_P(10,2,3).n(75)
2273            -7.194963600000000000000e8
2274        """
2275
2276        try:
2277            step_parent = kwds['parent']
2278        except KeyError:
2279            step_parent = parent(args[-1])
2280
2281        import sage.libs.mpmath.all as mpmath
2282
2283        try:
2284            precision = step_parent.prec()
2285        except AttributeError:
2286            precision = mpmath.mp.prec
2287
2288        return mpmath.call(
2289            mpmath.legenp,args[0],args[1],args[-1],prec = precision)
2290
2291    def _eval_special_values_(self,*args):
2292        """
2293        Special values known.
2294        EXAMPLES:
2295           
2296            sage: n, m = var('n m')
2297            sage: gen_legendre_P(n,m,0)
2298            2^m*gamma(1/2*m + 1/2*n + 1/2)*cos(1/2*(m + n)*pi)/(sqrt(pi)*gamma(-1/2*m + 1/2*n + 1))
2299        """
2300           
2301        if args[1] == 0:
2302            return legendre_P(args[0],args[-1])
2303
2304        if (args[-1] == 0):
2305            if is_Expression(args[0]): 
2306                try:
2307                    return cos(SR.pi()/2*(args[0]+args[1]))/(sqrt(SR.pi()))*\
2308                           gamma((args[0]+args[1]+1)/2)/\
2309                           gamma((args[0]-args[1])/2 + 1)*2**(args[1])
2310                except TypeError:
2311                    pass
2312            else:
2313                return cos(math.pi/2*(args[0]+args[1]))/(sqrt(math.pi))*\
2314                       gamma((args[0]+args[1]+1)/2)/\
2315                       gamma((args[0]-args[1])/2. + 1)*2**args[1]
2316
2317        raise ValueError("Value not found")
2318
2319    def _maxima_init_evaled_(self, *args):
2320        n = args[0]
2321        m = args[1]
2322        x = args[2]
2323        if is_Expression(n) or is_Expression(m):
2324            return None
2325       
2326        from sage.functions.all import sqrt
2327   
2328        if m.mod(2).is_zero() or m.is_one():
2329            return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'\
2330                                         %(ZZ(n),ZZ(m))), locals={'x':x})
2331        else:
2332            return sqrt(1-x**2)*(((n-m+1)*x*gen_legendre_P(n,m-1,x)-\
2333                                  (n+m-1)*gen_legendre_P(n-1,m-1,x))/(1-x**2))
2334
2335
2336    def _eval_numpy_(self, *args):
2337        #TODO: numpy_eval with help of the new scipy version!!!!
2338        #Reason scipy 0.8 supports stable and fast numerical evaluation
2339        #of ortho polys
2340        #Now this only evaluates the array pointwise, and only the first one...
2341        #if isinstance(arg[0], numpy.ndarray).
2342        #This is a hack to provide compability!
2343        """
2344        EXAMPLES::
2345        sage: import numpy
2346        sage: z = numpy.array([1,2])
2347        sage: z2 = numpy.array([[1,2],[1,2]])
2348        sage: z3 = numpy.array([1,2,3.])
2349        sage: gen_legendre_P(1,1,z/5.)
2350        array([-0.9797959 , -0.91651514])
2351        sage: gen_legendre_P(1,1,z2/5.)
2352        array([[-0.9797959 , -0.91651514],
2353               [-0.9797959 , -0.91651514]])
2354        sage: gen_legendre_P(1,1,z3/5.)
2355        array([-0.9797959 , -0.91651514, -0.8       ])
2356        """
2357
2358        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
2359            raise NotImplementedError(
2360           "Support for numpy array in first argument(s) is not supported yet!")
2361       
2362        result = numpy.zeros(args[-1].shape).tolist()
2363        if isinstance(args[-1][0],int):
2364            for k in xrange(len(args[-1])):
2365                result[k] = gen_legendre_P(args[0],args[1],ZZ(args[-1][k]))
2366
2367        if isinstance(args[-1][0],float):
2368            for k in xrange(len(args[-1])):
2369                result[k] = gen_legendre_P(args[0],args[1],RR(args[-1][k]))
2370
2371        if isinstance(args[-1][0],numpy.ndarray):
2372            for k in xrange(len(args[-1])):
2373                result[k] = gen_legendre_P(args[0],args[1],args[-1][k])
2374
2375        return numpy.array(result)
2376     
2377    def _derivative_(self,*args,**kwds):
2378        """return the derivative of gen_legendre_P in
2379        form of the Legendre Function.
2380        EXAMPLES::
2381        sage: n,m = var('n m')
2382        sage: derivative(gen_legendre_P(n,m,x),x)
2383        (n*x*gen_legendre_P(n, m, x) - (m + n)*gen_legendre_P(n - 1, m, x))/(x^2 - 1)
2384        sage: derivative(gen_legendre_P(3,1,x),x)
2385        -15*sqrt(-x^2 + 1)*x + 3/2*(5*(x - 1)^2 + 10*x - 6)*x/sqrt(-x^2 + 1)
2386        """
2387        diff_param = kwds['diff_param']
2388        if diff_param in [0,1]: 
2389            raise NotImplementedError(
2390            "Derivative w.r.t. to the index is not supported,\
2391            yet, and perhaps never will be...")
2392        else:
2393            return (args[0]*args[-1]*gen_legendre_P(args[0],args[1],args[-1])\
2394                    -(args[0]+args[1])*\
2395                    gen_legendre_P(args[0]-1,args[1],args[-1]))/(args[-1]**2-1)
2396
2397gen_legendre_P = Func_gen_legendre_P()
2398
2399class Func_gen_legendre_Q(OrthogonalPolynomial):
2400
2401    def __init__(self):
2402        OrthogonalPolynomial.__init__(self,"gen_legendre_Q",nargs = 3,
2403conversions =dict(maxima='assoc_legendre_q',mathematica='LegendreQ'))
2404
2405    def _evalf_(self, *args,**kwds):
2406        """
2407        Evals gen_legendre_Q
2408        numerically with mpmath.
2409        EXAMPLES::
2410            sage: gen_legendre_Q(10,2,3).n(75)
2411            -2.773909528741569374688e-7 - 1.130182239430298584113e9*I
2412        """
2413
2414        try:
2415            step_parent = kwds['parent']
2416        except KeyError:
2417            step_parent = parent(args[-1])
2418
2419        import sage.libs.mpmath.all as mpmath
2420
2421        try:
2422            precision = step_parent.prec()
2423        except AttributeError:
2424            precision = mpmath.mp.prec
2425
2426        return mpmath.call(
2427            mpmath.legenq,args[0],args[1],args[-1],prec = precision)
2428
2429    def _eval_special_values_(self,*args):
2430        """
2431        Special values known.
2432        EXAMPLES:
2433           
2434            sage: n, m = var('n m')
2435            sage: gen_legendre_Q(n,m,0)
2436            -sqrt(pi)*2^(m - 1)*gamma(1/2*m + 1/2*n + 1/2)*sin(1/2*(m + n)*pi)/gamma(-1/2*m + 1/2*n + 1)
2437        """
2438           
2439        if args[1] == 0:
2440            return legendre_Q(args[0],args[-1])
2441
2442        if (args[-1] == 0):
2443            if is_Expression(args[0]): 
2444                try:
2445                    return -(sqrt(SR.pi()))*sin(SR.pi()/2*(args[0]+args[1]))\
2446                           *gamma((args[0]+args[1]+1)/2)/\
2447                           gamma((args[0]-args[1])/2 + 1)*2**(args[1]-1)
2448                except TypeError:
2449                    pass
2450            else:
2451                return -(sqrt(math.pi))/2*sin(math.pi/2*(args[0]+args[1]))\
2452                       *gamma((args[0]+args[1]+1)/2)/\
2453                       gamma((args[0]-args[1])/2. + 1)*2**args[1]
2454
2455        raise ValueError("Value not found")
2456
2457    def _maxima_init_evaled_(self, *args):
2458        n = args[0]
2459        m = args[1]
2460        x = args[2]
2461        if is_Expression(n) or is_Expression(m):
2462            return None
2463       
2464        from sage.functions.all import sqrt
2465   
2466        if m <= n:
2467            return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'\
2468                                         %(ZZ(n),ZZ(m))), locals={'x':x})
2469        if m == n + 1 or n == 0:
2470            if m.mod(2).is_zero():
2471                denom = (1 - x**2)**(m/2)
2472            else:
2473                denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)   
2474                if m == n + 1:
2475                    return (-1)**m*(m-1).factorial()*2**n/denom
2476                else:
2477                    return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
2478        else:
2479            return ((n-m+1)*x*gen_legendre_Q(n,m-1,x)-\
2480                    (n+m-1)*gen_legendre_Q(n-1,m-1,x))/sqrt(1-x**2)
2481
2482
2483    def _eval_numpy_(self, *args):
2484        #TODO: numpy_eval with help of the new scipy version!!!!
2485        #Reason scipy 0.8 supports stable and fast numerical evaluation
2486        #of ortho polys
2487        #Now this only evaluates the array pointwise, and only the first one...
2488        #if isinstance(arg[0], numpy.ndarray).
2489        #This is a hack to provide compability!
2490        """
2491        EXAMPLES::
2492        sage: import numpy
2493        sage: z = numpy.array([1,2])
2494        sage: z2 = numpy.array([[1,2],[1,2]])
2495        sage: z3 = numpy.array([1,2,3.])
2496        sage: gen_legendre_Q(1,1,z/5.)
2497        array([-0.40276067, -0.82471644])
2498        sage: gen_legendre_Q(1,1,z2/5.)
2499        array([[-0.40276067, -0.82471644],
2500               [-0.40276067, -0.82471644]])
2501        sage: gen_legendre_Q(1,1,z3/5.)
2502        array([-0.40276067, -0.82471644, -1.30451774])
2503        """
2504
2505        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
2506            raise NotImplementedError(
2507           "Support for numpy array in first argument(s) is not supported yet!")
2508       
2509        result = numpy.zeros(args[-1].shape).tolist()
2510        if isinstance(args[-1][0],int):
2511            for k in xrange(len(args[-1])):
2512                result[k] = gen_legendre_Q(args[0],args[1],ZZ(args[-1][k]))
2513
2514        if isinstance(args[-1][0],float):
2515            for k in xrange(len(args[-1])):
2516                result[k] = gen_legendre_Q(args[0],args[1],RR(args[-1][k]))
2517
2518        if isinstance(args[-1][0],numpy.ndarray):
2519            for k in xrange(len(args[-1])):
2520                result[k] = gen_legendre_Q(args[0],args[1],args[-1][k])
2521
2522        return numpy.array(result)
2523     
2524    def _derivative_(self,*args,**kwds):
2525        """return the derivative of gen_legendre_Q in
2526        form of the Legendre Function.
2527        EXAMPLES::
2528        sage: n,m = var('n m')
2529        sage: derivative(gen_legendre_Q(n,m,x),x)
2530        (n*x*gen_legendre_Q(n, m, x) - (m + n)*gen_legendre_Q(n - 1, m, x))/(x^2 - 1)
2531        sage: derivative(gen_legendre_Q(0,1,x),x)
2532        -x/(-x^2 + 1)^(3/2)
2533        """
2534        diff_param = kwds['diff_param']
2535        if diff_param in [0,1]: 
2536            raise NotImplementedError("Derivative w.r.t. to the index \
2537            is not supported, yet, and perhaps never will be...")
2538        else:
2539            return (args[0]*args[-1]*gen_legendre_Q(args[0],args[1],args[-1])\
2540                    -(args[0]+args[1])*\
2541                    gen_legendre_Q(args[0]-1,args[1],args[-1]))/(args[-1]**2-1)
2542
2543
2544gen_legendre_Q = Func_gen_legendre_Q()
2545
2546
2547
2548
2549
2550
2551
2552
2553
2554
2555
2556