# Ticket #9706: orthogonal_polys.py

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

A new version of the orthogonal_polys.py file.

Line
1r"""
2Orthogonal Polynomials
3
4This module wraps some of the orthogonal/special functions in the
5Maxima package "orthopoly". This package was written by Barton
6Willis of the University of Nebraska at Kearney. It is released
7under the terms of the General Public License (GPL). Send
8Maxima-related bug reports and comments on this module to
10version information.
11
12
13-  The Chebyshev polynomial of the first kind arises as a solution
14   to the differential equation
15
16   .. math::
17
18         (1-x^2)\,y'' - x\,y' + n^2\,y = 0
19
20
21   and those of the second kind as a solution to
22
23   .. math::
24
25         (1-x^2)\,y'' - 3x\,y' + n(n+2)\,y = 0.
26
27
28   The Chebyshev polynomials of the first kind are defined by the
29   recurrence relation
30
31   .. math::
32
33     T_0(x) = 1 \, T_1(x) = x \, T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x). \,
34
35
36   The Chebyshev polynomials of the second kind are defined by the
37   recurrence relation
38
39   .. math::
40
41     U_0(x) = 1 \, U_1(x) = 2x \, U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x). \,
42
43
44
45   For integers m,n, they satisfy the orthogonality
46   relations
47
48   .. math::
49
50     \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.
51
52
53   and
54
55
56   .. math::
57
58     \int_{-1}^1 U_n(x)U_m(x)\sqrt{1-x^2}\,dx =\frac{\pi}{2}\delta_{m,n}.
59
60
61
62   They are named after Pafnuty Chebyshev (alternative
63   transliterations: Tchebyshef or Tschebyscheff).
64
65-  The Hermite polynomials are defined either by
66
67   .. math::
68
69     H_n(x)=(-1)^n e^{x^2/2}\frac{d^n}{dx^n}e^{-x^2/2}
70
71
72   (the "probabilists' Hermite polynomials"), or by
73
74
75   .. math::
76
77     H_n(x)=(-1)^n e^{x^2}\frac{d^n}{dx^n}e^{-x^2}
78
79
80   (the "physicists' Hermite polynomials"). Sage (via Maxima)
81   implements the latter flavor. These satisfy the orthogonality
82   relation
83
84   .. math::
85
86     \int_{-\infty}^\infty H_n(x)H_m(x)\,e^{-x^2}\,dx ={n!2^n}{\sqrt{\pi}}\delta_{nm}
87
88
89
90   They are named in honor of Charles Hermite.
91
92-  Each *Legendre polynomial* P_n(x) is an n-th degree polynomial.
93   It may be expressed using Rodrigues' formula:
94
95   .. math::
96
97      P_n(x) = (2^n n!)^{-1} {\frac{d^n}{dx^n} } \left[ (x^2 -1)^n \right].
98
99   These are solutions to Legendre's differential equation:
100
101   .. math::
102
103      {\frac{d}{dx}} \left[ (1-x^2) {\frac{d}{dx}} P(x) \right] + n(n+1)P(x) = 0.
104
105   and satisfy the orthogonality relation
106
107   .. math::
108
109      \int_{-1}^{1} P_m(x) P_n(x)\,dx = {\frac{2}{2n + 1}} \delta_{mn}
110
111   The *Legendre function of the second kind* Q_n(x) is another
112   (linearly independent) solution to the Legendre differential equation.
113   It is not an "orthogonal polynomial" however.
114
115   The associated Legendre functions of the first kind
116   P_\ell^m(x) can be given in terms of the "usual"
117   Legendre polynomials by
118
119   .. math::
120
121     \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}
122
123
124   Assuming 0 \le m \le \ell, they satisfy the orthogonality
125   relation:
126
127   .. math::
128
129      \int_{-1}^{1} P_k ^{(m)} P_\ell ^{(m)} dx  = \frac{2 (\ell+m)!}{(2\ell+1)(\ell-m)!}\ \delta _{k,\ell},
130
131
132   where \delta _{k,\ell} is the Kronecker delta.
133
134   The associated Legendre functions of the second kind
135   Q_\ell^m(x) can be given in terms of the "usual"
136   Legendre polynomials by
137
138
139   .. math::
140
141     Q_\ell^m(x)   =  (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}Q_\ell(x).
142
143
144
145   They are named after Adrien-Marie Legendre.
146
147-  Laguerre polynomials may be defined by the Rodrigues formula
148
149   .. math::
150
151      L_n(x)=\frac{e^x}{n!}\frac{d^n}{dx^n}\left(e^{-x} x^n\right).
152
153
154   They are solutions of Laguerre's equation:
155
156
157   .. math::
158
159      x\,y'' + (1 - x)\,y' + n\,y = 0\,
160
161   and satisfy the orthogonality relation
162
163
164   .. math::
165
166      \int_0^\infty L_m(x) L_n(x) e^{-x}\,dx = \delta_{mn}.
167
168
169
170   The generalized Laguerre polynomials may be defined by the
171   Rodrigues formula:
172
173
174   .. math::
175
176       L_n^{(\alpha)}(x)   = {\frac{x^{-\alpha} e^x}{n!}}{\frac{d^n}{dx^n}} \left(e^{-x} x^{n+\alpha}\right) .
177
178
179   (These are also sometimes called the associated Laguerre
180   polynomials.) The simple Laguerre polynomials are recovered from
181   the generalized polynomials by setting \alpha =0.
182
183   They are named after Edmond Laguerre.
184
185-  Jacobi polynomials are a class of orthogonal polynomials. They
186   are obtained from hypergeometric series in cases where the series
187   is in fact finite:
188
189   .. math::
190
191     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) ,
192
193
194   where ()_n is Pochhammer's symbol (for the rising
195   factorial), (Abramowitz and Stegun p561.) and thus have the
196   explicit expression
197
198
199   .. math::
200
201     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 .
202
203
204
205   They are named after Carl Jacobi.
206
207-  Ultraspherical or Gegenbauer polynomials are given in terms of
208   the Jacobi polynomials P_n^{(\alpha,\beta)}(x) with
209   \alpha=\beta=a-1/2 by
210
211
212   .. math::
213
214     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).
215
216
217   They satisfy the orthogonality relation
218
219   .. math::
220
221     \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)} ,
222
223
224   for a>-1/2. They are obtained from hypergeometric series
225   in cases where the series is in fact finite:
226
227
228   .. math::
229
230     C_n^{(a)}(z) =\frac{(2a)^{\underline{n}}}{n!} \,_2F_1\left(-n,2a+n;a+\frac{1}{2};\frac{1-z}{2}\right)
231
232
233   where \underline{n} is the falling factorial. (See
234   Abramowitz and Stegun p561)
235
236   They are named for Leopold Gegenbauer (1849-1903).
237
238
239For completeness, the Pochhammer symbol, introduced by Leo August
240Pochhammer, (x)_n, is used in the theory of special
241functions to represent the "rising factorial" or "upper factorial"
242
243.. math::
244
245         (x)_n=x(x+1)(x+2)\cdots(x+n-1)=\frac{(x+n-1)!}{(x-1)!}.
246
247
248On the other hand, the "falling factorial" or "lower factorial" is
249
250.. math::
251
252     x^{\underline{n}}=\frac{x!}{(x-n)!} ,
253
254
255in the notation of Ronald L. Graham, Donald E. Knuth and Oren
256Patashnik in their book Concrete Mathematics.
257
258.. note::
259
260   The first call of any of these will usually cost a bit extra
261   (it loads "specfun", but I'm not sure if that is the real reason).
262   The next call is usually faster but not always.
263
264TODO: Implement associated Legendre polynomials and Zernike
265polynomials. (Neither is in Maxima.)
266http://en.wikipedia.org/wiki/Associated_Legendre_polynomials
267http://en.wikipedia.org/wiki/Zernike_polynomials
268
269REFERENCES:
270
271-  Abramowitz and Stegun: Handbook of Mathematical Functions,
272   http://www.math.sfu.ca/ cbm/aands/
273
274-  http://en.wikipedia.org/wiki/Chebyshev_polynomials
275
276-  http://en.wikipedia.org/wiki/Legendre_polynomials
277
278-  http://en.wikipedia.org/wiki/Hermite_polynomials
279
280-  http://mathworld.wolfram.com/GegenbauerPolynomial.html
281
282-  http://en.wikipedia.org/wiki/Jacobi_polynomials
283
284-  http://en.wikipedia.org/wiki/Laguerre_polynomia
285
286-  http://en.wikipedia.org/wiki/Associated_Legendre_polynomials
287
288
289AUTHORS:
290
291- David Joyner (2006-06)
292- Stefan Reiterer (2010-)
293"""
294
295#*****************************************************************************
296#       Copyright (C) 2006 William Stein <wstein@gmail.com>
297#                     2006 David Joyner <wdj@usna.edu>
298#
300#
301#    This code is distributed in the hope that it will be useful,
302#    but WITHOUT ANY WARRANTY; without even the implied warranty of
303#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
304#    General Public License for more details.
305#
306#  The full text of the GPL is available at:
307#
309#*****************************************************************************
310
311import copy
312import sage.plot.plot
313import sage.interfaces.all
314from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
315from sage.rings.rational_field import RationalField
316from sage.rings.real_mpfr import RealField
317from sage.misc.sage_eval import sage_eval
318from sage.rings.all import QQ, ZZ, CDF, RDF
319import sage.rings.commutative_ring as commutative_ring
320import sage.rings.ring as ring
321from sage.calculus.calculus import maxima
322from sage.symbolic.function import BuiltinFunction, GinacFunction
323from sage.symbolic.expression import is_Expression
324import sage.functions.special
325from sage.functions.special import MaximaFunction, meval
326import numpy,scipy
327import math
328import sage.libs.mpmath.all as mpmath
329
330_done = False
331def _init():
332    """
333    Internal function which checks if Maxima has loaded the
334    "orthopoly" package.  All functions using this in this
335    file should call this function first.
336
337    TEST:
338
339    The global starts False::
340
341        sage: sage.functions.orthogonal_polys._done
342        False
343
344    Then after using one of these functions, it changes::
345
346        sage: from sage.functions.orthogonal_polys import chebyshev_T
347        sage: chebyshev_T(2,x)
348        2*(x - 1)^2 + 4*x - 3
349        sage: sage.functions.orthogonal_polys._done
350        True
351
352    Note that because here we use a Pynac variable x,
353    the representation of the function is different from
354    its actual doctest, where a polynomial indeterminate
355    x is used.
356    """
357    global _done
358    if _done:
359        return
361    # TODO -- make it possible to use the intervals returned
363    maxima.eval('orthopoly_returns_intervals:false;')
364    _done = True
365
367
368class OrthogonalPolynomial(BuiltinFunction):
369    """
370    Base Class for Orthogonal Polynomials. The evaluation as a polynomial
371    is done via maxima due performance reasons. Therefore the internal name
372    in maxima maxima_name has to be declared.
373    Convention: The first argument is always the order of the polynomial,
374    he last one is always the value x where the polynomial is evaluated.
375
376    """
377    def __init__(self, name, maxima_name, nargs = 2, latex_name = None, conversions = None):
378        self._maxima_name = maxima_name
379        BuiltinFunction.__init__(self, name = name, nargs = nargs, latex_name = latex_name, conversions = conversions)
380
381    def _maxima_init_evaled_(self, *args):
382        """
383        Returns a string which represents this function evaluated at
384        *args* in Maxima.
385        Remark: This function is stolen from the class MaximaFunction,
386        from sage.functions.special.py
387
388        EXAMPLES::
389
390         sage: chebyshev_T(3,x)
391         sage: 4*(x - 1)^3 + 12*(x - 1)^2 + 9*x - 8
392        """
393        args_maxima = []
394        for a in args:
395            if isinstance(a, str):
396                args_maxima.append(a)
397            elif hasattr(a, '_maxima_init_'):
398                args_maxima.append(a._maxima_init_())
399            else:
400                args_maxima.append(str(a))
401        return "%s(%s)"%(self._maxima_name, ', '.join(args_maxima))
402
403    def _clenshaw_method_(self,*args):
404        """
405        The Clenshaw method uses the three term recursion of the polynomial,
406        or explicit formulas instead of maxima to evaluate the polynomial
407        efficiently, if the x argument is not a symbolic expression.
408        The name comes from the Clenshaw algorithm for fast evaluation of
409        polynomialsin chebyshev form.
410
411        comparison: Maxima and Clenshaw algorithm for non-symbolic evaluation:
412         sage: time chebyshev_T(50,10) #clenshaw
413         CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
414         Wall time: 0.00 s
415         49656746733678312490954442369580252421769338391329426325400124999
416         sage: time sage.functions.orthogonal_polys.chebyshev_T(50,10)
417         #maxima
418         CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
419         Wall time: 0.05 s
420         49656746733678312490954442369580252421769338391329426325400124999
421
422         sage: time chebyshev_T(500,10); #clenshaw
423         CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
424         Wall time: 0.00 s
425         sage: time sage.functions.orthogonal_polys.chebyshev_T(500,10);
426         #maxima
427         CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s
428         Wall time: 0.77 s
429        """
430        raise NotImplementedError(
431            "No recursive calculation of values implemented (yet)!")
432
433    def _eval_special_values_(self,*args):
434        """
435        Evals the polynomial explicitly for special values.
436        EXAMPLES:
437        sage: chebyshev_T(n,-1)
438        (-1)^n
439        """
440        raise ValueError("No special values known!")
441
442    def _eval_(self, *args):
443        """
444
445        The symbolic evaluation is done with maxima, because the evaluation of
446        the Polynomial representation seems to be quite clever.
447        For the fast numerical evaluation an other method should be used...
448        Therefore I suggest Clenshaw's algorithm, which uses the rekursion!
449
450        EXAMPLES::
451
452        sage: time chebyshev_T(5,x)
453        CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
454        Wall time: 0.16 s
455        16*(x - 1)^5 + 80*(x - 1)^4 + 140*(x - 1)^3 + 100*(x - 1)^2 + 25*x - 24
456
457        time chebyshev_T(50,x)
458        CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
459        Wall time: 0.20 s
460        562949953421312*(x - 1)^50 + 28147497671065600*(x - 1)^49 +....
461
462         time chebyshev_T(100,x);
463         CPU times: user 0.12 s, sys: 0.00 s, total: 0.12 s
464         Wall time: 0.30 s
465
466         The function also checks for special values, and if
467         the order is an integer and in range!
468         EXAMPLES::
469         sage: var('n')
470         n
471         sage: chebyshev_T(n,-1)
472         (-1)^n
473         sage: chebyshev_T(-7,x)
474         ...
475         TypeError: Order is not valid!
476         sage: chebyshev_T(7.5,x)
477         TypeError: Order is not an Integer!
478
479        """
480        if not is_Expression(args[0]):
481            if args[0] != floor(args[0]):
482                raise TypeError("Order is not an Integer!")
483
484            if args[0] < -1:
485                raise TypeError("Order is not valid!")
486
487        try:
488            return self._eval_special_values_(*args)
489        except ValueError:
490            pass
491
492        if not (is_Expression(args[0]) or is_Expression(args[-1])):
493            try:
494                 return self._clenshaw_method_(*args)
495            except NotImplementedError:
496                pass
497
498        _init()
499        try:
500            s = maxima(self._maxima_init_evaled_(*args))
501        except TypeError:
502            return None
503        if self._maxima_name in repr(s):
504            return None
505        else:
506            return s.sage()
507
508
509        #TODO: Numerical Evals with help of the new scipy version!!!!
510        #Reason scipy 0.8 supports stable and fast numerical evaluation
511        #of ortho polys
512
513
514class Func_chebyshev_T(OrthogonalPolynomial):
515
516    """
517    Class for the Chebyshev polynomial of the first kind for integers
518    n>-1.
519
520    REFERENCE:
521
522    - AS 22.5.31 page 778 and AS 6.1.22 page 256.
523
524    EXAMPLES::
525
526       sage: chebyshev_T(3,x)
527       sage: 4*(x - 1)^3 + 12*(x - 1)^2 + 9*x - 8
528       sage: var('k')
529       k
530       sage: test = chebyshev_T(k,x)
531       sage: test
532       chebyshev_T(k, x)
533    """
534
535    def __init__(self):
536        OrthogonalPolynomial.__init__(self,"chebyshev_T","chebyshev_t",nargs = 2,
537conversions =dict(maxima='chebyshev_t',mathematica='ChebyshevT'))
538
539    def _eval_special_values_(self,*args):
540        """
541        Values known at the boundary.
542        EXAMPLES:
543        sage: chebyshev_T(n,1)
544        1
545        sage: chebyshev_T(n,-1)
546        (-1)^n
547        """
548        if args[-1] == 1:
549            return 1
550
551        if args[-1] == -1:
552            return (-1)**args[0]
553
555
556    def _clenshaw_method_(self,*args):
557            """
558            Clenshaw method for chebyshev_T (means use recursions in this case)
559            This is much faster for numerical evaluation than maxima!
560            """
561
562            k = args[0]
563            x = args[1]
564
565            if k == -1:
566                return 0
567            elif k == 0:
568                return 1
569            elif k == 1:
570                return x
571            else:
572                #TODO: When evaluation of Symbolic Expressions works better
573                #use these explicit formulas instead!
574                #if -1 <= x <= 1:
575                #    return cos(k*acos(x))
576                #elif 1 < x:
577                #    return cosh(k*acosh(x))
578                #else: # x < -1
579                #    return (-1)**(k%2)*cosh(k*acosh(-x))
580
581                help1 = 1
582                help2 = x
583                for j in xrange(0,k-1):
584                    help3 = 2*x*help2 - help1
585                    help1 = help2
586                    help2 = help3
587
588                return help3
589
590    def _evalf_(self, *args,**kwds):
591        """
592        Evals chebyshev_T numerically with mpmath.
593        EXAMPLES::
594        chebyshev_T(10,3).n(75)
595        2.261953700000000000000e7
596        """
597        parent = kwds['parent']
598        return mpmath.call(mpmath.chebyt,args[0],args[-1],prec = parent.prec())
599
600    def _derivative_(self, *args, **kwds):
601        """
602        Returns the derivative of chebyshev_T in form of the chebyshev Polynomial
603        of the second kind chebyshev_U
604        EXAMPLES::
605            sage: derivative(chebyshev_T(k,x),x)
606            k*chebyshev_U(k - 1, x)
607            sage: derivative(chebyshev_T(3,x),x)
608            12*(x - 1)^2 + 24*x - 15
609            sage: derivative(chebyshev_T(k,x),k)
610            ...
611            NotImplementedError: Derivative w.r.t. to the index is not
612            supported, yet, and perhaps never will be...
613
614        """
615        diff_param = kwds['diff_param']
616        if diff_param == 0: #args[0] is args[diff_param]:
617            raise NotImplementedError(
618"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
619        else:
620            return args[0]*chebyshev_U(args[0]-1,args[1])
621
622chebyshev_T = Func_chebyshev_T()
623
624class Func_chebyshev_U(OrthogonalPolynomial):
625
626    """
627    Class for the Chebyshev polynomial of the second kind for integers n>=-1.
628
629    REFERENCE:
630
631    - AS, 22.8.3 page 783 and AS 6.1.22 page 256.
632
633    EXAMPLES::
634
635        sage: chebyshev_U(3,x)
636        8*(x - 1)^3 + 24*(x - 1)^2 + 20*x - 16
637    """
638
639    def __init__(self):
640        OrthogonalPolynomial.__init__(self,"chebyshev_U","chebyshev_u",
641nargs = 2,conversions =dict(maxima='chebyshev_u',mathematica='ChebyshevU'))
642
643    def _clenshaw_method_(self,*args):
644            """
645            Clenshaw method for chebyshev_T (means use the recursion...)
646            This is much faster for numerical evaluation than maxima!
647            """
648
649            k = args[0]
650            x = args[1]
651
652            if k == -1:
653                return 0
654            elif k == 0:
655                return 1
656            elif k == 1:
657                return 2*x
658            else:
659                help1 = 1
660                help2 = 2*x
661                for j in xrange(0,k-1):
662                    help3 = 2*x*help2 - help1
663                    help1 = help2
664                    help2 = help3
665
666                return help3
667
668    def _eval_special_values_(self,*args):
669        """
670        Values known at the boundary.
671        EXAMPLES:
672        sage: chebyshev_U(n,1)
673        n + 1
674        sage: chebyshev_U(n,-1)
675        (n + 1)*(-1)^n
676        """
677        if args[-1] == 1:
678            return (args[0]+1)
679
680        if args[-1] == -1:
681            return (-1)**args[0]*(args[0]+1)
682
684
685    def _evalf_(self, *args,**kwds):
686        """
687        Evals chebyshev_U numerically with mpmath.
688        EXAMPLES::
689        chebyshev_U(10,3).n(75)
690        4.661117900000000000000e7
691        """
692        parent = kwds['parent']
693        return mpmath.call(mpmath.chebyu,args[0],args[-1],prec = parent.prec())
694
695    def _derivative_(self, *args, **kwds):
696        """
697        Returns the derivative of chebyshev_U in form of the chebyshev
698        Polynomials of the first and second kind
699        EXAMPLES::
700            sage: derivative(chebyshev_U(k,x),x)
701            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
702            sage: derivative(chebyshev_U(3,x),x)
703            24*(x - 1)^2 + 48*x - 28
704            sage: derivative(chebyshev_U(k,x),k)
705            ...
706            NotImplementedError: Derivative w.r.t. to the index is not
707            supported, yet, and perhaps never will be...
708
709        """
710        diff_param = kwds['diff_param']
711        if diff_param == 0: #args[0] is args[diff_param]:
712            raise NotImplementedError(
713"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
714        else:
715            return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
716                    chebyshev_U(args[0],args[1]))/(args[1]**2-1)
717
718chebyshev_U = Func_chebyshev_U()
719
720def gen_laguerre(n,a,x):
721    """
722    Returns the generalized Laguerre polynomial for integers n > -1.
723    Typically, a = 1/2 or a = -1/2.
724
725    REFERENCE:
726
727    - table on page 789 in AS.
728
729    EXAMPLES::
730
731        sage: x = PolynomialRing(QQ, 'x').gen()
732        sage: gen_laguerre(2,1,x)
733        1/2*x^2 - 3*x + 3
734        sage: gen_laguerre(2,1/2,x)
735        1/2*x^2 - 5/2*x + 15/8
736        sage: gen_laguerre(2,-1/2,x)
737        1/2*x^2 - 3/2*x + 3/8
738        sage: gen_laguerre(2,0,x)
739        1/2*x^2 - 2*x + 1
740        sage: gen_laguerre(3,0,x)
741        -1/6*x^3 + 3/2*x^2 - 3*x + 1
742    """
743    from sage.functions.all import sqrt
744    _init()
745    return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
746
747def gen_legendre_P(n,m,x):
748    r"""
749    Returns the generalized (or associated) Legendre function of the
750    first kind for integers n > -1, m > -1.
751
752    The awkward code for when m is odd and 1 results from the fact that
753    Maxima is happy with, for example, (1 - t^2)^3/2, but
754    Sage is not. For these cases the function is computed from the
755    (m-1)-case using one of the recursions satisfied by the Legendre
756    functions.
757
758    REFERENCE:
759
760    - Gradshteyn and Ryzhik 8.706 page 1000.
761
762    EXAMPLES::
763
764        sage: P.<t> = QQ[]
765        sage: gen_legendre_P(2, 0, t)
766        3/2*t^2 - 1/2
767        sage: gen_legendre_P(2, 0, t) == legendre_P(2, t)
768        True
769        sage: gen_legendre_P(3, 1, t)
770        -3/2*sqrt(-t^2 + 1)*(5*t^2 - 1)
771        sage: gen_legendre_P(4, 3, t)
772        105*sqrt(-t^2 + 1)*(t^2 - 1)*t
773        sage: gen_legendre_P(7, 3, I).expand()
774        -16695*sqrt(2)
775        sage: gen_legendre_P(4, 1, 2.5)
776        -583.562373654533*I
777    """
778    from sage.functions.all import sqrt
779    _init()
780    if m.mod(2).is_zero() or m.is_one():
781        return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
782    else:
783        return sqrt(1-x**2)*(((n-m+1)*x*gen_legendre_P(n,m-1,x)-(n+m-1)*gen_legendre_P(n-1,m-1,x))/(1-x**2))
784
785def gen_legendre_Q(n,m,x):
786    """
787    Returns the generalized (or associated) Legendre function of the
788    second kind for integers n>-1, m>-1.
789
790    Maxima restricts m = n. Hence the cases m n are computed using the
791    same recursion used for gen_legendre_P(n,m,x) when m is odd and
792    1.
793
794    EXAMPLES::
795
796        sage: P.<t> = QQ[]
797        sage: gen_legendre_Q(2,0,t)
798        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
799        sage: gen_legendre_Q(2,0,t) - legendre_Q(2, t)
800        0
801        sage: gen_legendre_Q(3,1,0.5)
802        2.49185259170895
803        sage: gen_legendre_Q(0, 1, x)
804        -1/sqrt(-x^2 + 1)
805        sage: gen_legendre_Q(2, 4, x).factor()
806        48*x/((x - 1)^2*(x + 1)^2)
807    """
808    from sage.functions.all import sqrt
809    if m <= n:
810        _init()
811        return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
812    if m == n + 1 or n == 0:
813        if m.mod(2).is_zero():
814            denom = (1 - x**2)**(m/2)
815        else:
816            denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)
817        if m == n + 1:
818            return (-1)**m*(m-1).factorial()*2**n/denom
819        else:
820            return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
821    else:
822        return ((n-m+1)*x*gen_legendre_Q(n,m-1,x)-(n+m-1)*gen_legendre_Q(n-1,m-1,x))/sqrt(1-x**2)
823
824def hermite(n,x):
825    """
826    Returns the Hermite polynomial for integers n > -1.
827
828    REFERENCE:
829
830    - AS 22.5.40 and 22.5.41, page 779.
831
832    EXAMPLES::
833
834        sage: x = PolynomialRing(QQ, 'x').gen()
835        sage: hermite(2,x)
836        4*x^2 - 2
837        sage: hermite(3,x)
838        8*x^3 - 12*x
839        sage: hermite(3,2)
840        40
841        sage: S.<y> = PolynomialRing(RR)
842        sage: hermite(3,y)
843        8.00000000000000*y^3 - 12.0000000000000*y
844        sage: R.<x,y> = QQ[]
845        sage: hermite(3,y^2)
846        8*y^6 - 12*y^2
847        sage: w = var('w')
848        sage: hermite(3,2*w)
849        8*(8*w^2 - 3)*w
850    """
851    _init()
852    return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
853
854def jacobi_P(n,a,b,x):
855    r"""
856    Returns the Jacobi polynomial P_n^{(a,b)}(x) for
857    integers n > -1 and a and b symbolic or a > -1
858    and b > -1. The Jacobi polynomials are actually defined
859    for all a and b. However, the Jacobi polynomial weight
860    (1-x)^a(1+x)^b isn't integrable for a \leq -1
861    or b \leq -1.
862
863    REFERENCE:
864
865    - table on page 789 in AS.
866
867    EXAMPLES::
868
869        sage: x = PolynomialRing(QQ, 'x').gen()
870        sage: jacobi_P(2,0,0,x)
871        3/2*x^2 - 1/2
872        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
873        5.009999999999998
874    """
875    _init()
876    return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
877
878def laguerre(n,x):
879    """
880    Returns the Laguerre polynomial for integers n > -1.
881
882    REFERENCE:
883
884    - AS 22.5.16, page 778 and AS page 789.
885
886    EXAMPLES::
887
888        sage: x = PolynomialRing(QQ, 'x').gen()
889        sage: laguerre(2,x)
890        1/2*x^2 - 2*x + 1
891        sage: laguerre(3,x)
892        -1/6*x^3 + 3/2*x^2 - 3*x + 1
893        sage: laguerre(2,2)
894        -1
895    """
896    _init()
897    return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
898
899def legendre_P(n,x):
900    """
901    Returns the Legendre polynomial of the first kind for integers
902    n > -1.
903
904    REFERENCE:
905
906    - AS 22.5.35 page 779.
907
908    EXAMPLES::
909
910        sage: P.<t> = QQ[]
911        sage: legendre_P(2,t)
912        3/2*t^2 - 1/2
913        sage: legendre_P(3, 1.1)
914        1.67750000000000
915        sage: legendre_P(3, 1 + I)
916        7/2*I - 13/2
917        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
918        [-179  242]
919        [-484  547]
920        sage: legendre_P(3, GF(11)(5))
921        8
922    """
923    _init()
924    return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)), locals={'x':x})
925
926def legendre_Q(n,x):
927    """
928    Returns the Legendre function of the second kind for integers
929    n>-1.
930
931    Computed using Maxima.
932
933    EXAMPLES::
934
935        sage: P.<t> = QQ[]
936        sage: legendre_Q(2, t)
937        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
938        sage: legendre_Q(3, 0.5)
939        -0.198654771479482
940        sage: legendre_Q(4, 2)
941        443/16*I*pi + 443/16*log(3) - 365/12
942        sage: legendre_Q(4, 2.0)
943        0.00116107583162324 + 86.9828465962674*I
944    """
945    _init()
946    return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
947
948def ultraspherical(n,a,x):
949    """
950    Returns the ultraspherical (or Gegenbauer) polynomial for integers
951    n > -1.
952
953    Computed using Maxima.
954
955    REFERENCE:
956
957    - AS 22.5.27
958
959    EXAMPLES::
960
961        sage: x = PolynomialRing(QQ, 'x').gen()
962        sage: ultraspherical(2,3/2,x)
963        15/2*x^2 - 3/2
964        sage: ultraspherical(2,1/2,x)
965        3/2*x^2 - 1/2
966        sage: ultraspherical(1,1,x)
967        2*x
968        sage: t = PolynomialRing(RationalField(),"t").gen()
969        sage: gegenbauer(3,2,t)
970        32*t^3 - 12*t
971    """
972    _init()
973    return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
974
975gegenbauer = ultraspherical