# Ticket #9706: orthogonal_polys.5.py

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

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