# Ticket #9706: orthogonal_polys.4.py

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

Latest version. It holds classes of all polys (but not all completed yet)

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, is_inexact
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
329from sage.functions.other import floor, gamma, factorial, abs, binomial
330from sage.functions.trig import sin
331from sage.functions.log import ln
332import sage.symbolic.expression as expression
333from sage.structure.parent import Parent
334from sage.structure.coerce import parent
335
336_done = False
337def _init():
338    """
339    Internal function which checks if Maxima has loaded the
340    "orthopoly" package.  All functions using this in this
341    file should call this function first.
342
343    TEST:
344
345    The global starts False::
346
347        sage: sage.functions.orthogonal_polys._done
348        False
349
350    Then after using one of these functions, it changes:: (The value is now
351    False for chebyshev_T because chebyshev_T uses clenshaw method instead...)
352
353        sage: from sage.functions.orthogonal_polys import chebyshev_T
354        sage: chebyshev_T(2,x)
355        2*x^2 - 1
356        sage: sage.functions.orthogonal_polys._done
357        False
358        sage: legendre_Q(1,0.1)
359        -0.989966465226892
360        sage: sage.functions.orthogonal_polys._done
361        True
362
363    Note that because here we use a Pynac variable x,
364    the representation of the function is different from
365    its actual doctest, where a polynomial indeterminate
366    x is used.
367    """
368    global _done
369    if _done:
370        return
372    # TODO -- make it possible to use the intervals returned
374    maxima.eval('orthopoly_returns_intervals:false;')
375    _done = True
376
378
379class OrthogonalPolynomial(BuiltinFunction):
380    """
381    Base Class for Orthogonal Polynomials. The evaluation as a polynomial
382    is done via maxima due performance reasons. Therefore the internal name
383    in maxima maxima_name has to be declared.
384    Convention: The first argument is always the order of the polynomial,
385    he last one is always the value x where the polynomial is evaluated.
386
387    """
388    def __init__(self, name, nargs = 2, latex_name = None, conversions = {}):
389        try:
390            self._maxima_name = conversions['maxima']
391        except KeyError:
392            self._maxima_name = None
393
394        BuiltinFunction.__init__(self, name = name, nargs = nargs, latex_name = latex_name, conversions = conversions)
395
396    def _maxima_init_evaled_(self, *args):
397        """
398        Returns a string which represents this function evaluated at
399        *args* in Maxima.
400        In fact these are thought to be the old wrappers for the orthogonal
401        polynomials. These are used when the other evaluation methods fail,
402        or are not fast enough. Experiments showed that for the symbolic
403        evaluation for larger n maxima is faster, but for small n simply use
404        of the recursion formulas is faster. A little switch does the trick...
405
406        EXAMPLES::
407
408            sage: chebyshev_T(3,x)
409            4*x^3 - 3*x
410        """
411        return None
412
413    def _clenshaw_method_(self,*args):
414        """
415        The Clenshaw method uses the three term recursion of the polynomial,
416        or explicit formulas instead of maxima to evaluate the polynomial
417        efficiently, if the x argument is not a symbolic expression.
418        The name comes from the Clenshaw algorithm for fast evaluation of
419        polynomialsin chebyshev form.
420
421        comparison: Maxima and Clenshaw algorithm for non-symbolic evaluation:
422         #sage: time chebyshev_T(50,10) #clenshaw
423         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
424         #Wall time: 0.00 s
425         #49656746733678312490954442369580252421769338391329426325400124999
426         #sage: time sage.functions.orthogonal_polys.chebyshev_T(50,10)
427         #maxima
428         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
429         #Wall time: 0.05 s
430         #49656746733678312490954442369580252421769338391329426325400124999
431
432         #sage: time chebyshev_T(500,10); #clenshaw
433         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
434         #Wall time: 0.00 s
435         #sage: time sage.functions.orthogonal_polys.chebyshev_T(500,10);
436         #maxima
437         #CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s
438         #Wall time: 0.77 s
439        """
440        raise NotImplementedError(
441            "No recursive calculation of values implemented (yet)!")
442
443    def _eval_special_values_(self,*args):
444        """
445        Evals the polynomial explicitly for special values.
446        EXAMPLES:
447
448            sage: var('n')
449            n
450            sage: chebyshev_T(n,-1)
451            (-1)^n
452        """
453        raise ValueError("No special values known!")
454
455
456    def _eval_(self, *args):
457        """
458
459        The symbolic evaluation is done with maxima, because the evaluation of
460        the Polynomial representation seems to be quite clever.
461        For the fast numerical evaluation an other method should be used...
462        Therefore I suggest Clenshaw's algorithm, which uses the rekursion!
463        The function also checks for special values, and if
464        the order is an integer and in range!
465
466        Update: There are cases where the pure Clanshaw algorithm seems faster
467        and sometimes maxima seems faster. e.g. for legendre_P clenshaw is faster
468        for legendre_Q maxima is faster.
469
470        performance:
471            #sage: time chebyshev_T(5,x) #maxima
472            #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
473            #Wall time: 0.16 s
474            #16*(x - 1)^5 + 80*(x - 1)^4 + 140*(x - 1)^3 + 100*(x - 1)^2 + 25*x - 24
475
476            #sage: time chebyshev_T(5,x) #clenshaw
477            #CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
478            #Wall time: 0.01 s
479            #16*x^5 - 20*x^3 + 5*x
480
481            #time chebyshev_T(50,x)
482            #CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
483            #Wall time: 0.20 s
484            #562949953421312*(x - 1)^50 + 28147497671065600*(x - 1)^49 +....
485
486            #time chebyshev_T(100,x);
487            #CPU times: user 0.12 s, sys: 0.00 s, total: 0.12 s
488            #Wall time: 0.30 s
489
490        EXAMPLES::
491            sage: chebyshev_T(5,x)
492            16*x^5 - 20*x^3 + 5*x
493            sage: var('n')
494            n
495            sage: chebyshev_T(n,-1)
496            (-1)^n
497            sage: chebyshev_T(-7,x)
498            chebyshev_T(-7, x)
499            sage: chebyshev_T(3/2,x)
500            chebyshev_T(3/2, x)
501
502        """
503
504        if not is_Expression(args[0]):
505
506            if not is_Expression(args[-1]) and is_inexact(args[-1]):
507                try:
508                    return self._evalf_(*args)
509                except AttributeError:
510                    pass
511                except mpmath.NoConvergence:
512                    print "Warning: mpmath returns NoConvergence!"
513                    print "Switching to clenshaw_method, but it may not be stable!"
514
515
516            #A faster check would be nice...
517            if args[0] != floor(args[0]):
518                if not is_Expression(args[-1]):
519                    try:
520                        return self._evalf_(*args)
521                    except AttributeError:
522                        pass
523                else:
524                    return None
525
526            if args[0] < 0:
527                return None
528
529
530        try:
531            return self._eval_special_values_(*args)
532        except ValueError:
533            pass
534
535
536        if not is_Expression(args[0]):
537
538            try:
539                return self._clenshaw_method_(*args)
540            except NotImplementedError:
541                pass
542
543        if self._maxima_name is None:
544            return None
545        else:
546            _init()
547            try:
548                #s = maxima(self._maxima_init_evaled_(*args))
549                #This above is very inefficient! The older
550                #methods were much faster...
551                return self._maxima_init_evaled_(*args)
552            except TypeError:
553                return None
554            if self._maxima_name in repr(s):
555                return None
556            else:
557                return s.sage()
558
559
560        #TODO: numpy_eval with help of the new scipy version!!!!
561        #Reason scipy 0.8 supports stable and fast numerical evaluation
562        #of ortho polys
563
564
565class Func_chebyshev_T(OrthogonalPolynomial):
566
567    """
568    Class for the Chebyshev polynomial of the first kind.
569
570    REFERENCE:
571
572    - AS 22.5.31 page 778 and AS 6.1.22 page 256.
573
574    EXAMPLES::
575
576       sage: chebyshev_T(3,x)
577       4*x^3 - 3*x
578       sage: var('k')
579       k
580       sage: test = chebyshev_T(k,x)
581       sage: test
582       chebyshev_T(k, x)
583    """
584
585    def __init__(self):
586        OrthogonalPolynomial.__init__(self,"chebyshev_T",nargs = 2,
587conversions =dict(maxima='chebyshev_t',mathematica='ChebyshevT'))
588
589    def _eval_special_values_(self,*args):
590        """
591        Values known for special values of x.
592        For details see A.S. 22.4 (p. 777)
593
594        EXAMPLES:
595
596            sage: var('n')
597            n
598            sage: chebyshev_T(n,1)
599            1
600            sage: chebyshev_T(n,-1)
601            (-1)^n
602        """
603        if args[-1] == 1:
604            return 1
605
606        if args[-1] == -1:
607            return (-1)**args[0]
608
609        if (args[-1] == 0):
610            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
611
613
614    def _evalf_(self, *args,**kwds):
615        """
616        Evals chebyshev_T
617        numerically with mpmath.
618        EXAMPLES::
619            sage: chebyshev_T(10,3).n(75)
620            2.261953700000000000000e7
621        """
622        try:
623            step_parent = kwds['parent']
624        except KeyError:
625            step_parent = parent(args[-1])
626
627        try:
628            precision = step_parent.prec()
629        except AttributeError:
630            precision = mpmath.mp.prec
631
632        return mpmath.call(mpmath.chebyt,args[0],args[-1],prec = precision)
633
634    def _maxima_init_evaled_(self, *args):
635        n = args[0]
636        x = args[1]
637        return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
638
639    def _clenshaw_method_(self,*args):
640            """
641            Clenshaw method for chebyshev_T (means use recursions in this case)
642            This is much faster for numerical evaluation than maxima!
643            See A.S. 227 (p. 782) for details for the recurions
644            """
645
646            k = args[0]
647            x = args[1]
648
649            if k == 0:
650                return 1
651            elif k == 1:
652                return x
653            else:
654                #TODO: When evaluation of Symbolic Expressions works better
655                #use these explicit formulas instead!
656                #if -1 <= x <= 1:
657                #    return cos(k*acos(x))
658                #elif 1 < x:
659                #    return cosh(k*acosh(x))
660                #else: # x < -1
661                #    return (-1)**(k%2)*cosh(k*acosh(-x))
662
663                help1 = 1
664                help2 = x
665                if is_Expression(x):
666                    #raise NotImplementedError
667                    help3 = 0
668                    for j in xrange(0,floor(k/2)+1):
669                        help3 = help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j-1)/factorial(j)/factorial(k-2*j)
670                    help3 = help3*k/2
671                else:
672                     for j in xrange(0,k-1):
673                        help3 = 2*x*help2 - help1
674                        help1 = help2
675                        help2 = help3
676
677                return help3
678
679    def _derivative_(self, *args, **kwds):
680        """
681        Returns the derivative of chebyshev_T in form of the chebyshev Polynomial
682        of the second kind chebyshev_U
683        EXAMPLES::
684            sage: var('k')
685            k
686            sage: derivative(chebyshev_T(k,x),x)
687            k*chebyshev_U(k - 1, x)
688            sage: derivative(chebyshev_T(3,x),x)
689            12*x^2 - 3
690            sage: derivative(chebyshev_T(k,x),k)
691            Traceback (most recent call last):
692            ...
693            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
694
695        """
696        diff_param = kwds['diff_param']
697        if diff_param == 0:
698            raise NotImplementedError(
699"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
700        else:
701            return args[0]*chebyshev_U(args[0]-1,args[1])
702
703chebyshev_T = Func_chebyshev_T()
704
705class Func_chebyshev_U(OrthogonalPolynomial):
706
707    """
708    Class for the Chebyshev polynomial of the second kind.
709
710    REFERENCE:
711
712    - AS, 22.8.3 page 783 and AS 6.1.22 page 256.
713
714    EXAMPLES::
715
716        sage: chebyshev_U(3,x)
717        8*x^3 - 4*x
718    """
719
720    def __init__(self):
721        OrthogonalPolynomial.__init__(self,"chebyshev_U",
722nargs = 2,conversions =dict(maxima='chebyshev_u',mathematica='ChebyshevU'))
723
724    def _clenshaw_method_(self,*args):
725        """
726        Clenshaw method for chebyshev_U (means use the recursion...)
727        This is much faster for numerical evaluation than maxima!
728        See A.S. 227 (p. 782) for details for the recurions
729        """
730        k = args[0]
731        x = args[1]
732
733        if k == 0:
734            return 1
735        elif k == 1:
736            return 2*x
737        else:
738            help1 = 1
739            help2 = 2*x
740            if is_Expression(x):
741                #raise NotImplementedError("Maxima is faster here!")
742                help3 = 0
743                for j in xrange(0,floor(k/2)+1):
744                    help3 = help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j)/factorial(j)/factorial(k-2*j)
745                #help3 = help3*k/2
746            else:
747                for j in xrange(0,k-1):
748                    help3 = 2*x*help2 - help1
749                    help1 = help2
750                    help2 = help3
751
752            return help3
753
754    def _maxima_init_evaled_(self, *args):
755        """
756        Uses
757        EXAMPLES::
758        sage: x = PolynomialRing(QQ, 'x').gen()
759        sage: chebyshev_T(2,x)
760        2*x^2 - 1
761        """
762        n = args[0]
763        x = args[1]
764        return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
765
766
767    def _evalf_(self, *args,**kwds):
768        """
769        Evals chebyshev_U
770        numerically with mpmath.
771        EXAMPLES::
772            sage: chebyshev_U(10,3).n(75)
773            4.661117900000000000000e7
774        """
775        try:
776            step_parent = kwds['parent']
777        except KeyError:
778            step_parent = parent(args[-1])
779
780        try:
781            precision = step_parent.prec()
782        except AttributeError:
783            precision = mpmath.mp.prec
784
785        return mpmath.call(mpmath.chebyu,args[0],args[-1],prec = precision)
786
787    def _eval_special_values_(self,*args):
788        """
789        Special values known. A.S. 22.4 (p.777).
790        EXAMPLES:
791
792            sage: var('n')
793            n
794            sage: chebyshev_U(n,1)
795            n + 1
796            sage: chebyshev_U(n,-1)
797            (n + 1)*(-1)^n
798        """
799        if args[-1] == 1:
800            return (args[0]+1)
801
802        if args[-1] == -1:
803            return (-1)**args[0]*(args[0]+1)
804
805        if (args[-1] == 0):
806            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
807
809
810    def _derivative_(self, *args, **kwds):
811        """
812        Returns the derivative of chebyshev_U in form of the chebyshev
813        Polynomials of the first and second kind
814
815        EXAMPLES::
816            sage: var('k')
817            k
818            sage: derivative(chebyshev_U(k,x),x)
819            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
820            sage: derivative(chebyshev_U(3,x),x)
821            24*x^2 - 4
822            sage: derivative(chebyshev_U(k,x),k)
823            Traceback (most recent call last):
824            ...
825            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
826
827        """
828        diff_param = kwds['diff_param']
829        if diff_param == 0:
830            raise NotImplementedError(
831"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
832        else:
833            return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
834                    chebyshev_U(args[0],args[1]))/(args[1]**2-1)
835
836chebyshev_U = Func_chebyshev_U()
837
838class Func_gegenbauer(OrthogonalPolynomial):
839    """
840    Returns the ultraspherical (or Gegenbauer) polynomial.
841
842    REFERENCE:
843
844    - AS 22.5.27
845
846    EXAMPLES::
847
848        sage: x = PolynomialRing(QQ, 'x').gen()
849        sage: ultraspherical(2,3/2,x)
850        15/2*x^2 - 3/2
851        sage: ultraspherical(2,1/2,x)
852        3/2*x^2 - 1/2
853        sage: ultraspherical(1,1,x)
854        2*x
855        sage: t = PolynomialRing(RationalField(),"t").gen()
856        sage: gegenbauer(3,2,t)
857        32*t^3 - 12*t
858    """
859    def __init__(self):
860        OrthogonalPolynomial.__init__(self,"gegenbauer",nargs = 3,conversions =dict(maxima='ultraspherical',mathematica='GegenbauerC'))
861
862    def _clenshaw_method_(self,*args):
863        """
864        Clenshaw method for gegenbauer poly (means use the recursion...)
865        This is much faster for numerical evaluation than maxima!
866        See A.S. 227 (p. 782) for details for the recurions
867        """
868        k = args[0]
869        x = args[-1]
870        alpha = args[1]
871
872        if is_Expression(alpha) or abs(alpha) > numpy.finfo(float).eps:
873            alpha_zero = False
874        else:
875            alpha_zero = True
876
877        if k == 0:
878            return 1
879        elif k == 1:
880            if not alpha_zero:
881                return 2*x*alpha
882            else:
883                return 2*x # It seems that maxima evals this the wrong way! (see A.S. 22.4 (p.777))
884        else:
885            help1 = 1
886            if alpha_zero:
887                help2 = 2*x
888            else:
889                help2 = 2*x*alpha
890
891            help3 = 0
892            if is_Expression(x):
893                for j in xrange(0,floor(k/2)+1):
894                    help3 = help3 + (-1)**j*gamma(alpha+k-j)/factorial(j)/factorial(k-2*j)*(2*x)**(k-2*j)
895            else:
896                for j in xrange(1,k):
897                    help3 = 2*(j+alpha)*x*help2 - (j+2*alpha-1)*help1
898                    help3 = help3/(j+1)
899                    help1 = help2
900                    help2 = help3
901
902            return help3
903
904    def _maxima_init_evaled_(self, *args):
905        """
906        Uses
907        EXAMPLES::
908        sage: x = PolynomialRing(QQ, 'x').gen()
909        sage: ultraspherical(2,1/2,x)
910        3/2*x^2 - 1/2
911        """
912        n = args[0]
913        a = args[1]
914        x = args[2]
915        return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
916
917
918    def _evalf_(self, *args,**kwds):
919        """
920        Evals gegenbauer
921        numerically with mpmath.
922        EXAMPLES::
923            sage: gegenbauer(10,2,3.).n(54)
924            5.25360702000000e8
925        """
926        try:
927            step_parent = kwds['parent']
928        except KeyError:
929            step_parent = parent(args[-1])
930
931        try:
932            precision = step_parent.prec()
933        except AttributeError:
934            precision = mpmath.mp.prec
935
936        return mpmath.call(mpmath.gegenbauer,args[0],args[1],args[-1],prec = precision)
937
938    def _eval_special_values_(self,*args):
939        """
940        Special values known. A.S. 22.4 (p.777)
941        EXAMPLES:
942
943            sage: var('n a')
944            (n, a)
945            sage: gegenbauer(n,1/2,x)
946            legendre_P(n, x)
947            sage: gegenbauer(n,0,x)
948            1/2*n*chebyshev_T(n, x)
949            sage: gegenbauer(n,1,x)
950            chebyshev_U(n, x)
951            sage: gegenbauer(n,a,1)
952            binomial(2*a + n - 1, n)
953            sage: gegenbauer(n,a,-1)
954            (-1)^n*binomial(2*a + n - 1, n)
955            sage: gegenbauer(n,a,0)
956            1/2*((-1)^n + 1)*(-1)^(1/2*n)*gamma(a + 1/2*n)/(gamma(a)*gamma(1/2*n))
957        """
958        if args[1] == 0 and args[0] != 0:
959            return args[0]*chebyshev_T(args[0],args[-1])/2
960
961        if args[1] == 0.5:
962            return legendre_P(args[0],args[-1])
963
964        if args[1] == 1:
965            return chebyshev_U(args[0],args[-1])
966
967        if args[-1] == 1:
968            return binomial(args[0] + 2*args[1] - 1,args[0])
969
970        if args[-1] == -1:
971            return (-1)**args[0]*binomial(args[0] + 2*args[1] - 1,args[0])
972
973        if (args[-1] == 0):
974            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2*gamma(args[1]+args[0]/2)/gamma(args[1])/gamma(args[0]/2)
975
977
978    def _derivative_(self, *args, **kwds):
979        """
980        Returns the derivative of chebyshev_U in form of the chebyshev
981        Polynomials of the first and second kind
982
983        EXAMPLES::
984            sage: var('k a')
985            (k, a)
986            sage: derivative(gegenbauer(k,a,x),x)
987            (k*x*gegenbauer(k, a, x) - (2*a + k - 1)*gegenbauer(k - 1, a, x))/(x^2 - 1)
988            sage: derivative(gegenbauer(4,3,x),x)
989            1920*x^3 - 480*x
990            sage: derivative(gegenbauer(k,a,x),a)
991            Traceback (most recent call last):
992            ...
993            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
994
995        """
996        diff_param = kwds['diff_param']
997        if diff_param in [0,1]:
998            raise NotImplementedError(
999"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1000        else:
1001            return (-args[0]*args[-1]*gegenbauer(args[0],args[1],args[2])+ (args[0] + 2*args[1]-1)*gegenbauer(args[0]-1,args[1],args[2]))/(1-args[-1]**2)
1002
1003gegenbauer = Func_gegenbauer()
1004
1005class Func_gen_laguerre(OrthogonalPolynomial):
1006    """
1007    Returns the generalized Laguerre polynomial for integers n > -1.
1008    Typically, a = 1/2 or a = -1/2.
1009
1010    REFERENCE:
1011
1012    - table on page 789 in AS.
1013
1014    EXAMPLES::
1015
1016        sage: x = PolynomialRing(QQ, 'x').gen()
1017        sage: gen_laguerre(2,1,x)
1018        1/2*x^2 - 3*x + 3
1019        sage: gen_laguerre(2,1/2,x)
1020        1/2*x^2 - 5/2*x + 15/8
1021        sage: gen_laguerre(2,-1/2,x)
1022        1/2*x^2 - 3/2*x + 3/8
1023        sage: gen_laguerre(2,0,x)
1024        1/2*x^2 - 2*x + 1
1025        sage: gen_laguerre(3,0,x)
1026        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1027    """
1028    def __init__(self):
1029        OrthogonalPolynomial.__init__(self,"gen_laguerre",nargs = 3,
1030conversions =dict(maxima='gen_laguerre',mathematica='LaguerreL'))
1031
1032    def _clenshaw_method_(self,*args):
1033        """
1034        Clenshaw method for gen_laguerre polynomial (means use the recursion...)
1035        This is much faster for numerical evaluation than maxima!
1036        See A.S. 227 (p. 782) for details for the recurions.
1037        Warning: The clanshaw method for the laguerre Polynomials
1038        should only used for exact data types, when high orders are
1039        used, due to weak instabilities of the recursion!
1040        """
1041        k = args[0]
1042        a = args[1]
1043        x = args[-1]
1044
1045        if k == 0:
1046            return 1
1047        elif k == 1:
1048            return 1-x + a
1049        else:
1050            help1 = 1
1051            help2 = 1-x + a
1052            if is_Expression(x):
1053                help3 = 0
1054                for j in xrange(0,k+1):
1055                    help3 = help3 + (-1)**j*binomial(k+a,k-j)/factorial(j)*x**j
1056            else:
1057                for j in xrange(1,k):
1058                    help3 = -x*help2 + (2*j+1+a)*help2 - (j+a)*help1
1059                    help3 = help3/(j+1)
1060                    help1 = help2
1061                    help2 = help3
1062
1063            return help3
1064
1065    def _maxima_init_evaled_(self, *args):
1066        n = args[0]
1067        a = args[1]
1068        x = args[-1]
1069
1070        from sage.functions.all import sqrt
1071        _init()
1072        return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
1073
1074    def _evalf_(self, *args,**kwds):
1075        """
1076        Evals laguerre polynomial
1077        numerically with mpmath.
1078        EXAMPLES::
1079            sage: laguerre(3,5.).n(53)
1080            2.66666666666667
1081        """
1082        try:
1083            step_parent = kwds['parent']
1084        except KeyError:
1085            step_parent = parent(args[-1])
1086
1087        try:
1088            precision = step_parent.prec()
1089        except AttributeError:
1090            precision = mpmath.mp.prec
1091
1092        return mpmath.call(mpmath.laguerre,args[0],args[1],args[-1],prec = precision)
1093
1094    def _eval_special_values_(self,*args):
1095        """
1096        Special values known.
1097        EXAMPLES:
1098
1099            sage: var('n')
1100            n
1101            sage: laguerre(n,0)
1102            1
1103        """
1104
1105        if (args[-1] == 0):
1106            try:
1107                return binomial(args[0]+args[1],args[0])
1108            except TypeError:
1109                pass
1110
1112
1113
1114    def _derivative_(self,*args,**kwds):
1115        """return the derivative of laguerre in
1116        form of the Laguerre Polynomial.
1117        EXAMPLES::
1118        sage: n = var('n')
1119        sage: derivative(laguerre(3,x),x)
1120        -1/2*x^2 + 3*x - 3
1121        sage: derivative(laguerre(n,x),x)
1122        -(n*laguerre(n - 1, x) - n*laguerre(n, x))/x
1123        """
1124        diff_param = kwds['diff_param']
1125        if diff_param == 0:
1126            raise NotImplementedError(
1127"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1128        else:
1129            return (args[0]*laguerre(args[0],args[1],args[-1])-(args[0]+args[1])*laguerre(args[0]-1,args[1],args[-1]))/args[-1]
1130
1131gen_laguerre = Func_gen_laguerre()
1132
1133class Func_hermite(OrthogonalPolynomial):
1134    """
1135    Returns the Hermite polynomial for integers n > -1.
1136
1137    REFERENCE:
1138
1139    - AS 22.5.40 and 22.5.41, page 779.
1140
1141    #sage: S.<y> = PolynomialRing(RR)<- Here a bug of the base class BuiltinFunction occours!
1142
1143    EXAMPLES::
1144
1145        sage: x = PolynomialRing(QQ, 'x').gen()
1146        sage: hermite(2,x)
1147        4*x^2 - 2
1148        sage: hermite(3,x)
1149        8*x^3 - 12*x
1150        sage: hermite(3,2)
1151        40
1152        sage: y = var('y')
1153        sage: hermite(3,y).polynomial(PolynomialRing(RR, 'y'))
1154        8.00000000000000*y^3 - 12.0000000000000*y
1155        sage: R.<x,y> = QQ[]
1156        sage: hermite(3,y^2)
1157        8*y^6 - 12*y^2
1158        sage: w = var('w')
1159        sage: hermite(3,2*w)
1160        64*w^3 - 24*w
1161    """
1162
1163    def __init__(self):
1164        OrthogonalPolynomial.__init__(self,"hermite",nargs = 2,
1165conversions =dict(maxima='hermite',mathematica='HermiteH'))
1166
1167    def _clenshaw_method_(self,*args):
1168        """
1169        Clenshaw method for hermite polynomial (means use the recursion...)
1170        See A.S. 227 (p. 782) for details for the recurions
1171        For the symbolic evaluation, maxima seems to be quite fast.
1172        The break even point between the recursion and Maxima is about
1173        n = 25
1174        """
1175        k = args[0]
1176        x = args[1]
1177
1178        if k == 0:
1179            return 1
1180        elif k == 1:
1181            return 2*x
1182        else:
1183            help1 = 1
1184            help2 = 2*x
1185            if is_Expression(x):
1186                help3 = 0
1187                for j in xrange(0,floor(k/2)+1):
1188                    help3 = help3 + (-1)**j*(2*x)**(k-2*j)/factorial(j)/factorial(k-2*j)
1189                help3 = help3*factorial(k)
1190            else:
1191                for j in xrange(1,k):
1192                    help3 = 2*x*help2 - 2*j*help1
1193                    help1 = help2
1194                    help2 = help3
1195
1196            return help3
1197
1198
1199    def _evalf_(self, *args,**kwds):
1200        """
1201        Evals hermite
1202        numerically with mpmath.
1203        EXAMPLES::
1204            sage: hermite(3,2.).n(74)
1205            40.0000000000000000000
1206        """
1207        try:
1208            step_parent = kwds['parent']
1209        except KeyError:
1210            step_parent = parent(args[-1])
1211
1212        try:
1213            precision = step_parent.prec()
1214        except AttributeError:
1215            precision = mpmath.mp.prec
1216
1217        return mpmath.call(mpmath.hermite,args[0],args[-1],prec = precision)
1218
1219    def _eval_special_values_(self,*args):
1220        """
1221        Special values known. A.S. 22.4 (p.777)
1222        EXAMPLES:
1223
1224            sage: var('n')
1225            n
1226            sage: hermite(n,0)
1227            ((-1)^n + 1)*(-1)^(1/2*n)*factorial(n)/gamma(1/2*n + 1)
1228        """
1229
1230        if (args[-1] == 0):
1231            return (1+(-1)**args[0])*(-1)**(args[0]/2)*factorial(args[0])/gamma(args[0]/2+1)
1232
1234
1235    def _maxima_init_evaled_(self, *args):
1236        """
1237        Old maxima method.
1238        """
1239        n = args[0]
1240        x = args[1]
1241        return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
1242
1243    def _derivative_(self, *args, **kwds):
1244        """
1245        Returns the derivative of the hermite polynomial in form of the chebyshev
1246        Polynomials of the first and second kind
1247
1248        EXAMPLES::
1249            sage: var('k')
1250            k
1251            sage: derivative(hermite(k,x),x)
1252            2*k*hermite(k - 1, x)
1253
1254        """
1255        diff_param = kwds['diff_param']
1256        if diff_param == 0:
1257            raise NotImplementedError(
1258"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1259        else:
1260            return 2*args[0]*hermite(args[0]-1,args[1])
1261
1262hermite = Func_hermite()
1263
1264
1265class Func_jacobi_P(OrthogonalPolynomial):
1266    r"""
1267    Returns the Jacobi polynomial P_n^{(a,b)}(x) for
1268    integers n > -1 and a and b symbolic or a > -1
1269    and b > -1. The Jacobi polynomials are actually defined
1270    for all a and b. However, the Jacobi polynomial weight
1271    (1-x)^a(1+x)^b isn't integrable for a \leq -1
1272    or b \leq -1.
1273
1274    REFERENCE:
1275
1276    - table on page 789 in AS.
1277
1278    EXAMPLES::
1279
1280        sage: x = PolynomialRing(QQ, 'x').gen()
1281        sage: jacobi_P(2,0,0,x)
1282        3/2*x^2 - 1/2
1283        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
1284        5.009999999999998
1285    """
1286
1287    def __init__(self):
1288        OrthogonalPolynomial.__init__(self,"jacobi_P",nargs = 4,
1289conversions =dict(maxima='jacobi_p',mathematica='JacobiP'))
1290
1291    def _clenshaw_method_(self,*args):
1292        """
1293        Clenshaw method for jacobi_P (means use the recursion,
1294        or sum formula)
1295        This is much faster for numerical evaluation than maxima!
1296        See A.S. 227 (p. 782) for details for the recurions.
1297        Warning: The clanshaw method for the Jacobi Polynomials
1298        should only used for exact data types, when high orders are
1299        used, due to weak instabilities of the recursion!
1300        """
1301        k = args[0]
1302        x = args[-1]
1303        alpha = args[1]
1304        beta = args[2]
1305
1306        if k == 0:
1307            return 1
1308        elif k == 1:
1309            return (alpha-beta + (alpha+beta+2)*x)/2
1310        else:
1311
1312            if is_Expression(x) or is_Expression(alpha) or is_Expression(beta):
1313                #Here we use the sum formula of jacobi_P it seems this is rather
1314                #optimal for use.
1315                help1 = gamma(alpha+k+1)/factorial(k)/gamma(alpha+beta+k+1)
1316                help2 = 0
1317                for j in xrange(0,k+1):
1318                    help2 = help2 + binomial(k,j)*gamma(alpha+beta+k+j+1)/gamma(alpha+j+1)*((x-1)/2)**j
1319                return help1*help2
1320            else:
1321                help1 = 1
1322                help2 = (alpha-beta + (alpha+beta+2)*x)/2
1323
1324                for j in xrange(1,k):
1325                    a1 = 2*(j+1)*(j+alpha+beta+1)*(2*j+alpha+beta)
1326                    a2 = (2*j+alpha+beta+1)*(alpha**2-beta**2)
1327                    a3 = (2*j+alpha+beta)*(2*j+alpha+beta+1)*(2*j+alpha+beta+2)
1328                    a4 = 2*(j+alpha)*(j+beta)*(2*j+alpha+beta+2)
1329                    help3 = (a2+a3*x)*help2 - a4*help1
1330                    help3 = help3/a1
1331                    help1 = help2
1332                    help2 = help3
1333
1334            return help3
1335
1336    def _maxima_init_evaled_(self, *args):
1337        """
1338        Old maxima method.
1339        """
1340        n = args[0]
1341        a = args[1]
1342        b = args[2]
1343        x = args[3]
1344        return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
1345
1346    def _evalf_(self, *args,**kwds):
1347        """
1348        Evals jacobi_P
1349        numerically with mpmath.
1350        EXAMPLES::
1351            sage: jacobi_P(10,2,3,3).n(75)
1352            1.322776620000000000000e8
1353        """
1354        try:
1355            step_parent = kwds['parent']
1356        except KeyError:
1357            step_parent = parent(args[-1])
1358
1359        try:
1360            precision = step_parent.prec()
1361        except AttributeError:
1362            precision = mpmath.mp.prec
1363
1364        return mpmath.call(mpmath.jacobi,args[0],args[1],args[2],args[-1],prec = precision)
1365
1366    def _eval_special_values_(self,*args):
1367        """
1368        Special values known. A.S. 22.4 (p.777)
1369        EXAMPLES:
1370
1371            sage: var('n')
1372            n
1373            sage: legendre_P(n,1)
1374            1
1375            sage: legendre_P(n,-1)
1376            (-1)^n
1377        """
1378        if args[-1] == 1:
1379            return binomial(args[0]+args[1],args[0])
1380
1381        if args[-1] == -1:
1382            return (-1)**args[0]*binomial(args[0]+args[1],args[0])
1383
1384        if args[1] == 0 and args[2] == 0:
1385            return legendre_P(args[0],args[-1])
1386
1387        if args[1] == -0.5 and args[2] == -0.5:
1388            try:
1389                return binomial(2*args[0],args[0])*chebyshev_T(args[0],args[-1])/4**args[0]
1390            except TypeError:
1391                pass
1392
1394
1395    def _derivative_(self, *args, **kwds):
1396        """
1397        Returns the derivative of jacobi_P in form of jacobi_polynomials
1398
1399        EXAMPLES::
1400            sage: var('k a b')
1401            (k, a, b)
1402            sage: derivative(jacobi_P(k,a,b,x),x)
1403            -(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))
1404            sage: derivative(jacobi_P(2,1,3,x),x)
1405            14*x - 7/2
1406            sage: derivative(jacobi_P(k,a,b,x),a)
1407            Traceback (most recent call last):
1408            ...
1409            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
1410
1411        """
1412        diff_param = kwds['diff_param']
1413        if diff_param in [0,1,2]:
1414            raise NotImplementedError(
1415"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1416        else:
1417            return (args[0]*(args[1]-args[2]-(2*args[0]+args[1]+args[2])*args[-1])*
1418                    jacobi_P(args[0],args[1],args[2],args[-1])+ 2*(args[0]+args[1])*
1419                    (args[0]+args[2])*jacobi_P(args[0]-1,args[1],args[2],args[-1]))/(2*args[0]+args[1]+args[2])/(1-args[-1]**2)
1420
1421
1422jacobi_P = Func_jacobi_P()
1423
1424class Func_laguerre(OrthogonalPolynomial):
1425    """
1426    Returns the Laguerre polynomial for integers n > -1.
1427
1428    REFERENCE:
1429
1430    - AS 22.5.16, page 778 and AS page 789.
1431
1432    EXAMPLES::
1433
1434        sage: x = PolynomialRing(QQ, 'x').gen()
1435        sage: laguerre(2,x)
1436        1/2*x^2 - 2*x + 1
1437        sage: laguerre(3,x)
1438        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1439        sage: laguerre(2,2)
1440        -1
1441    """
1442    def __init__(self):
1443        OrthogonalPolynomial.__init__(self,"laguerre",nargs = 2,
1444conversions =dict(maxima='laguerre',mathematica='LaguerreL'))
1445
1446    def _clenshaw_method_(self,*args):
1447        """
1448        Clenshaw method for laguerre polynomial (means use the recursion...)
1449        This is much faster for numerical evaluation than maxima!
1450        See A.S. 227 (p. 782) for details for the recurions.
1451        Warning: The clanshaw method for the laguerre Polynomials
1452        should only used for exact data types, when high orders are
1453        used, due to weak instabilities of the recursion!
1454        """
1455        k = args[0]
1456        x = args[-1]
1457
1458        if k == 0:
1459            return 1
1460        elif k == 1:
1461            return 1-x
1462        else:
1463            help1 = 1
1464            help2 = 1-x
1465            if is_Expression(x):
1466                help3 = 0
1467                for j in xrange(0,k+1):
1468                    help3 = help3 + (-1)**j*binomial(k,k-j)/factorial(j)*x**j
1469            else:
1470                for j in xrange(1,k):
1471                    help3 = -x*help2 + (2*j+1)*help2 - j*help1
1472                    help3 = help3/(j+1)
1473                    help1 = help2
1474                    help2 = help3
1475
1476            return help3
1477
1478    def _maxima_init_evaled_(self, *args):
1479        n = args[0]
1480        x = args[1]
1481        return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
1482
1483    def _evalf_(self, *args,**kwds):
1484        """
1485        Evals laguerre polynomial
1486        numerically with mpmath.
1487        EXAMPLES::
1488            sage: laguerre(3,5.).n(53)
1489            2.66666666666667
1490        """
1491        try:
1492            step_parent = kwds['parent']
1493        except KeyError:
1494            step_parent = parent(args[-1])
1495
1496        try:
1497            precision = step_parent.prec()
1498        except AttributeError:
1499            precision = mpmath.mp.prec
1500
1501        return mpmath.call(mpmath.laguerre,args[0],0,args[-1],prec = precision)
1502
1503    def _eval_special_values_(self,*args):
1504        """
1505        Special values known.
1506        EXAMPLES:
1507
1508            sage: var('n')
1509            n
1510            sage: laguerre(n,0)
1511            1
1512        """
1513
1514        if (args[-1] == 0):
1515            try:
1516                return 1
1517            except TypeError:
1518                pass
1519
1521
1522
1523    def _derivative_(self,*args,**kwds):
1524        """return the derivative of laguerre in
1525        form of the Laguerre Polynomial.
1526        EXAMPLES::
1527        sage: n = var('n')
1528        sage: derivative(laguerre(3,x),x)
1529        -1/2*x^2 + 3*x - 3
1530        sage: derivative(laguerre(n,x),x)
1531        -(n*laguerre(n - 1, x) - n*laguerre(n, x))/x
1532        """
1533        diff_param = kwds['diff_param']
1534        if diff_param == 0:
1535            raise NotImplementedError(
1536"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1537        else:
1538            return (args[0]*laguerre(args[0],args[-1])-args[0]*laguerre(args[0]-1,args[-1]))/args[1]
1539
1540laguerre = Func_laguerre()
1541
1542class Func_legendre_P(OrthogonalPolynomial):
1543    """
1544    Returns the Legendre polynomial of the first kind..
1545
1546    REFERENCE:
1547
1548    - AS 22.5.35 page 779.
1549
1550    EXAMPLES::
1551
1552        sage: P.<t> = QQ[]
1553        sage: legendre_P(2,t)
1554        3/2*t^2 - 1/2
1555        sage: legendre_P(3, 1.1)
1556        1.67750000000000
1557        sage: legendre_P(3, 1 + I)
1558        7/2*I - 13/2
1559        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
1560        Traceback (most recent call last):
1561        ...
1562        TypeError: cannot coerce arguments: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Integer Ring to Symbolic Ring
1563        sage: legendre_P(3, GF(11)(5))
1564        8
1565    """
1566    def __init__(self):
1567        OrthogonalPolynomial.__init__(self,"legendre_P",nargs = 2,
1568conversions =dict(maxima='legendre_p',mathematica='LegendreP'))
1569
1570    def _clenshaw_method_(self,*args):
1571        """
1572        Clenshaw method for legendre_P (means use the recursion...)
1573        This is much faster for numerical evaluation than maxima!
1574        See A.S. 227 (p. 782) for details for the recurions.
1575        Warning: The clanshaw method for the Legendre Polynomials
1576        should only used for exact data types, when high orders are
1577        used, due to weak instabilities of the recursion!
1578        """
1579        k = args[0]
1580        x = args[-1]
1581
1582        if k == 0:
1583            return 1
1584        elif k == 1:
1585            return x
1586        else:
1587            help1 = 1
1588            help2 = x
1589            if is_Expression(x):
1590                #raise NotImplementedError("Maxima is faster here...")
1591                help3 = 0
1592                for j in xrange(0,floor(k/2)+1):
1593                    help3 = help3 + (-1)**j*x**(k-2*j)*binomial(k,j)*binomial(2*k-2*j,k)
1594                help3 = help3/2**k
1595
1596            else:
1597                for j in xrange(1,k):
1598                    help3 = (2*j+1)*x*help2 - j*help1
1599                    help3 = help3/(j+1)
1600                    help1 = help2
1601                    help2 = help3
1602
1603            return help3
1604
1605    def _maxima_init_evaled_(self, *args):
1606        n = args[0]
1607        x = args[1]
1608        return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)), locals={'x':x})
1609
1610    def _evalf_(self, *args,**kwds):
1611        """
1612        Evals legendre_P
1613        numerically with mpmath.
1614        EXAMPLES::
1615            sage: legendre_P(10,3).n(75)
1616            8.097453000000000000000e6
1617        """
1618        try:
1619            step_parent = kwds['parent']
1620        except KeyError:
1621            step_parent = parent(args[-1])
1622
1623        try:
1624            precision = step_parent.prec()
1625        except AttributeError:
1626            precision = mpmath.mp.prec
1627
1628        return mpmath.call(mpmath.legendre,args[0],args[-1],prec = precision)
1629
1630
1631    def _eval_special_values_(self,*args):
1632        """
1633        Special values known.
1634        EXAMPLES:
1635
1636            sage: var('n')
1637            n
1638            sage: legendre_P(n,1)
1639            1
1640            sage: legendre_P(n,-1)
1641            (-1)^n
1642        """
1643        if args[-1] == 1:
1644            return 1
1645
1646        if args[-1] == -1:
1647            return (-1)**args[0]
1648
1649        if (args[-1] == 0):
1650            try:
1651                return (1+(-1)**args[0])/2*binomial(args[0],args[0]/2)/4**(args[0]/2)
1652            except TypeError:
1653                pass
1654
1656
1657    def _derivative_(self,*args,**kwds):
1658        """return the derivative of legendre_P in
1659        form of the Legendre Polynomial.
1660        EXAMPLES::
1661        derivative(legendre_P(n,x),x)
1662        -(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x - 1)^2
1663        derivative(legendre_P(3,x),x)
1664        15/2*x^2 - 3/2
1665        """
1666        diff_param = kwds['diff_param']
1667        if diff_param == 0:
1668            raise NotImplementedError(
1669"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1670        else:
1671            return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*legendre_P(args[0],args[-1]))/(1-args[-1])**2
1672
1673
1674legendre_P = Func_legendre_P()
1675
1676class Func_legendre_Q(OrthogonalPolynomial):
1677    """
1678    Return the chebyshev function of the second kind.
1679
1680    REFERENCE:
1681
1682    - A.S. 8 (p. 332)
1683
1684    EXAMPLES::
1685
1686        sage: P.<t> = QQ[]
1687        sage: legendre_Q(2, t)
1688        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
1689        sage: legendre_Q(3, 0.5)
1690        -0.198654771479482
1691        sage: legendre_Q(4, 2)
1692        443/16*I*pi + 443/16*log(3) - 365/12
1693        sage: legendre_Q(4, 2.0)
1694        0.00116107583162324 + 86.9828465962674*I
1695
1696    """
1697    def __init__(self):
1698        OrthogonalPolynomial.__init__(self,"legendre_Q",nargs = 2,
1699conversions =dict(maxima='legendre_q',mathematica='LegendreQ'))
1700
1701    def _maxima_init_evaled_(self, *args):
1702        """
1703        Maxima seems just fine for legendre Q. So we use it here!
1704        """
1705        n = args[0]
1706        x = args[1]
1707        return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
1708
1709    def _clenshaw_method_(self,*args):
1710        """
1711        Clenshaw method for legendre_q (means use the recursion...)
1712        This is much faster for numerical evaluation than maxima!
1713        See A.S. 8.5.3 (p. 334) for details for the recurions.
1714        Warning: The clanshaw method for the Legendre fUNCTIONS
1715        should only used for exact data types, when high orders are
1716        used, due to weak instabilities of the recursion!
1717        """
1718        raise NotImplementedError("Function not ready yet...")
1719
1720        k = args[0]
1721        x = args[-1]
1722
1723        if k == 0:
1724            return ln((1+x)/(1-x))/2
1725        elif k == 1:
1726            return x/2*ln((1+x)/(1-x))-1
1727        else:
1728            if is_Expression(x):
1729                raise NotImplementedError("Maxima works fine here!")
1730                #it seems that the old method just works fine here...
1731                #raise NotImplementedError("clenshaw does not work well...")
1732            else:
1733                help1 = ln((1+x)/(1-x))/2
1734                help2 = x/2*ln((1+x)/(1-x))-1
1735
1736                for j in xrange(1,k):
1737                    help3 = (2*j+1)*x*help2 - j*help1
1738                    help3 = help3/(j+1)
1739                    help1 = help2
1740                    help2 = help3
1741
1742            return help3
1743
1744    def _eval_special_values_(self,*args):
1745        """
1746        Special values known.
1747        EXAMPLES:
1748
1749            sage: var('n')
1750            n
1751            sage: legendre_P(n,1)
1752            1
1753            sage: legendre_P(n,-1)
1754            (-1)^n
1755        """
1756        if args[-1] == 1:
1757            return NaN
1758
1759        if args[-1] == -1:
1760            return NaN
1761
1762        if (args[-1] == 0):
1763            if is_Expression(args[0]):
1764                try:
1765                    return -(SR.pi()**(QQ(1)/QQ(2)))/2*sin(SR.pi()/2*args[0])*gamma((args[0]+1)/2)*gamma(args[0]/2 + 1)
1766                except TypeError:
1767                    pass
1768            else:
1769                return -(math.pi**(QQ(1)/QQ(2)))/2*sin(math.pi/2*args[0])*gamma((args[0]+1)/2)*gamma(args[0]/2. + 1)
1770
1772
1773    def _evalf_(self, *args,**kwds):
1774        """
1775        Evals legendre_Q
1776        numerically with mpmath.
1777        Function does not work yet...
1778        EXAMPLES::
1779            sage: legendre_P(10,3).n(75)
1780            8.097453000000000000000e6
1781        """
1782        raise AttributeError("mpmath function does not work correctly!")
1783
1784    def _derivative_(self,*args,**kwds):
1785        """return the derivative of legendre_Q in
1786        form of the Legendre Function.
1787        EXAMPLES::
1788        derivative(legendre_P(n,x),x)
1789        -(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x - 1)^2
1790        derivative(legendre_P(3,x),x)
1791        15/2*x^2 - 3/2
1792        """
1793        diff_param = kwds['diff_param']
1794        if diff_param == 0:
1795            raise NotImplementedError(
1796"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1797        else:
1798            return (gen_legendre_Q(args[0],1,args[-1]))/sqrt(1-args[-1]**2)
1799
1800legendre_Q = Func_legendre_Q()
1801
1802def legendre_Q_old(n,x):
1803    """
1804    Returns the Legendre function of the second kind for integers
1805    n>-1. (old version)
1806
1807    Computed using Maxima.
1808
1809    EXAMPLES::
1810
1811        sage: P.<t> = QQ[]
1812        sage: legendre_Q(2, t)
1813        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
1814        sage: legendre_Q(3, 0.5)
1815        -0.198654771479482
1816        sage: legendre_Q(4, 2)
1817        443/16*I*pi + 443/16*log(3) - 365/12
1818        sage: legendre_Q(4, 2.0)
1819        0.00116107583162324 + 86.9828465962674*I
1820    """
1821    _init()
1822    return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
1823
1824
1825
1826def gen_laguerre_old(n,a,x):
1827    """
1828    Returns the generalized Laguerre polynomial for integers n > -1.
1829    Typically, a = 1/2 or a = -1/2.
1830
1831    REFERENCE:
1832
1833    - table on page 789 in AS.
1834
1835    EXAMPLES::
1836
1837        sage: x = PolynomialRing(QQ, 'x').gen()
1838        sage: gen_laguerre(2,1,x)
1839        1/2*x^2 - 3*x + 3
1840        sage: gen_laguerre(2,1/2,x)
1841        1/2*x^2 - 5/2*x + 15/8
1842        sage: gen_laguerre(2,-1/2,x)
1843        1/2*x^2 - 3/2*x + 3/8
1844        sage: gen_laguerre(2,0,x)
1845        1/2*x^2 - 2*x + 1
1846        sage: gen_laguerre(3,0,x)
1847        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1848    """
1849    from sage.functions.all import sqrt
1850    _init()
1851    return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
1852
1853class Func_gen_legendre_P(OrthogonalPolynomial):
1854
1855    def __init__(self):
1856        OrthogonalPolynomial.__init__(self,"gen_legendre_P",nargs = 3,
1857conversions =dict(maxima='assoc_legendre_p',mathematica='LegendreP'))
1858
1859    def _maxima_init_evaled_(self, *args):
1860        n = args[0]
1861        m = args[1]
1862        x = args[2]
1863        if is_Expression(n) or is_Expression(m):
1864            return None
1865
1866        from sage.functions.all import sqrt
1867
1868        if m.mod(2).is_zero() or m.is_one():
1869            return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
1870        else:
1871            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))
1872
1873
1874    def _derivative_(self,*args,**kwds):
1875        """return the derivative of gen_legendre_P in
1876        form of the Legendre Function.
1877        EXAMPLES::
1878        derivative(legendre_P(n,x),x)
1879        -(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x - 1)^2
1880        derivative(legendre_P(3,x),x)
1881        15/2*x^2 - 3/2
1882        """
1883        diff_param = kwds['diff_param']
1884        if diff_param == 0:
1885            raise NotImplementedError("Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1886        else:
1887            return ((-args[1]*args[-1]*gen_legendre_P(args[0],args[1],args[-1]))+(args*[0]+args[1])*(args[0]-args[1]-1)*sqrt(args[0]^2-1)*gen_legendre_P(args[0],args[1]-1,args[-1]))/(args[0]**2-1)
1888
1889gen_legendre_P = Func_gen_legendre_P()
1890
1891class Func_gen_legendre_Q(OrthogonalPolynomial):
1892
1893    def __init__(self):
1894        OrthogonalPolynomial.__init__(self,"gen_legendre_Q",nargs = 3,
1895conversions =dict(maxima='assoc_legendre_q',mathematica='LegendreQ'))
1896
1897    def _maxima_init_evaled_(self, *args):
1898        n = args[0]
1899        m = args[1]
1900        x = args[2]
1901        if is_Expression(n) or is_Expression(m):
1902            return None
1903
1904        from sage.functions.all import sqrt
1905
1906        if m <= n:
1907            return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
1908        if m == n + 1 or n == 0:
1909            if m.mod(2).is_zero():
1910                denom = (1 - x**2)**(m/2)
1911            else:
1912                denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)
1913                if m == n + 1:
1914                    return (-1)**m*(m-1).factorial()*2**n/denom
1915                else:
1916                    return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
1917        else:
1918            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)
1919
1920
1921    #def _derivative_(self,*args,**kwds):
1922    #    """return the derivative of gen_legendre_P in
1923    #    form of the Legendre Function.
1924    #    EXAMPLES::
1925    #    derivative(legendre_P(n,x),x)
1926    #    -(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x - 1)^2
1927    #    derivative(legendre_P(3,x),x)
1928    #    15/2*x^2 - 3/2
1929    #    """
1930    #    diff_param = kwds['diff_param']
1931    #    if diff_param == 0:
1932    #        raise NotImplementedError("Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1933    #    else:
1934    #        return ((-args[1]*args[-1]*gen_legendre_P(args[0],args[1],args[-1]))+(args*[0]+args[1])*(args[0]-args[1]-1)*sqrt(args[0]^2-1)*gen_legendre_P(args[0],args[1]-1,args[-1]))/(args[0]**2-1)
1935
1936gen_legendre_Q = Func_gen_legendre_Q()
1937
1938def gen_legendre_P_old(n,m,x):
1939    r"""
1940    Returns the generalized (or associated) Legendre function of the
1941    first kind for integers n > -1, m > -1.
1942
1943    The awkward code for when m is odd and 1 results from the fact that
1944    Maxima is happy with, for example, (1 - t^2)^3/2, but
1945    Sage is not. For these cases the function is computed from the
1946    (m-1)-case using one of the recursions satisfied by the Legendre
1947    functions.
1948
1949    REFERENCE:
1950
1951    - Gradshteyn and Ryzhik 8.706 page 1000.
1952
1953    EXAMPLES::
1954
1955        sage: P.<t> = QQ[]
1956        sage: expand(gen_legendre_P(2, 0, t))
1957        3/2*t^2 - 1/2
1958        sage: expand(gen_legendre_P(2, 0, t)) - legendre_P(2, t)
1959        0
1960        sage: gen_legendre_P(3, 1, t)
1961        -3/2*sqrt(-t^2 + 1)*(5*(t - 1)^2 + 10*t - 6)
1962        sage: gen_legendre_P(4, 3, t)
1963        15*sqrt(-t^2 + 1)*((t^2 - 1)*(7*(t - 1)^2 + 14*t - 8)*t - 6*(t^2 - 1)*t)/(t^2 - 1)
1964        sage: gen_legendre_P(7, 3, I).expand()
1965        -16695*sqrt(2)
1966        sage: gen_legendre_P(4, 1, 2.5)
1967        -583.562373654533*I
1968    """
1969    from sage.functions.all import sqrt
1970    _init()
1971    if m.mod(2).is_zero() or m.is_one():
1972        return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
1973    else:
1974        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))
1975
1976def gen_legendre_Q_old(n,m,x):
1977    """
1978    Returns the generalized (or associated) Legendre function of the
1979    second kind for integers n>-1, m>-1.
1980
1981    Maxima restricts m = n. Hence the cases m n are computed using the
1982    same recursion used for gen_legendre_P(n,m,x) when m is odd and
1983    1.
1984
1985    EXAMPLES::
1986
1987        sage: P.<t> = QQ[]
1988        sage: gen_legendre_Q(2,0,t)
1989        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
1990        sage: gen_legendre_Q(2,0,t) - legendre_Q(2, t)
1991        0
1992        sage: gen_legendre_Q(3,1,0.5)
1993        2.49185259170895
1994        sage: gen_legendre_Q(0, 1, x)
1995        -1/sqrt(-x^2 + 1)
1996        sage: gen_legendre_Q(2, 4, x).factor()
1997        48*x/((x - 1)^2*(x + 1)^2)
1998    """
1999    from sage.functions.all import sqrt
2000    if m <= n:
2001        _init()
2002        return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
2003    if m == n + 1 or n == 0:
2004        if m.mod(2).is_zero():
2005            denom = (1 - x**2)**(m/2)
2006        else:
2007            denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)
2008        if m == n + 1:
2009            return (-1)**m*(m-1).factorial()*2**n/denom
2010        else:
2011            return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
2012    else:
2013        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)
2014
2015
2016def jacobi_P_old(n,a,b,x):
2017    r"""
2018    Returns the Jacobi polynomial P_n^{(a,b)}(x) for
2019    integers n > -1 and a and b symbolic or a > -1
2020    and b > -1. The Jacobi polynomials are actually defined
2021    for all a and b. However, the Jacobi polynomial weight
2022    (1-x)^a(1+x)^b isn't integrable for a \leq -1
2023    or b \leq -1.
2024
2025    REFERENCE:
2026
2027    - table on page 789 in AS.
2028
2029    EXAMPLES::
2030
2031        sage: x = PolynomialRing(QQ, 'x').gen()
2032        sage: jacobi_P(2,0,0,x)
2033        3/2*x^2 - 1/2
2034        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
2035        5.009999999999998
2036    """
2037    _init()
2038    return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
2039
2040def laguerre_old(n,x):
2041    """
2042    Returns the Laguerre polynomial for integers n > -1.
2043
2044    REFERENCE:
2045
2046    - AS 22.5.16, page 778 and AS page 789.
2047
2048    EXAMPLES::
2049
2050        sage: x = PolynomialRing(QQ, 'x').gen()
2051        sage: laguerre(2,x)
2052        1/2*x^2 - 2*x + 1
2053        sage: laguerre(3,x)
2054        -1/6*x^3 + 3/2*x^2 - 3*x + 1
2055        sage: laguerre(2,2)
2056        -1
2057    """
2058    _init()
2059    return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
2060
2061def ultraspherical(n,a,x):
2062    """
2063    Returns the ultraspherical (or Gegenbauer) polynomial for integers
2064    n > -1.
2065
2066    Computed using Maxima.
2067
2068    REFERENCE:
2069
2070    - AS 22.5.27
2071
2072    EXAMPLES::
2073
2074        sage: x = PolynomialRing(QQ, 'x').gen()
2075        sage: ultraspherical(2,3/2,x)
2076        15/2*x^2 - 3/2
2077        sage: ultraspherical(2,1/2,x)
2078        3/2*x^2 - 1/2
2079        sage: ultraspherical(1,1,x)
2080        2*x
2081        sage: t = PolynomialRing(RationalField(),"t").gen()
2082        sage: gegenbauer(3,2,t)
2083        32*t^3 - 12*t
2084    """
2085    _init()
2086    return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
2087
2088