Ticket #9706: orthogonal_polys.3.py

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

Version from 10. August 2010

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
9willisb@unk.edu. In your report, please include Maxima and specfun
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#
299#  Distributed under the terms of the GNU General Public License (GPL)
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#
308#                  http://www.gnu.org/licenses/
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
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: gegenbauer(1,2,3)
359        12
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
371    maxima.eval('load("orthopoly");')
372    # TODO -- make it possible to use the intervals returned
373    # instead of just discarding this info!
374    maxima.eval('orthopoly_returns_intervals:false;')
375    _done = True
376
377#load /home/maldun/sage/sage-4.5.1/devel/sage-ortho/sage/functions/orthogonal_polys.py
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            0
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
612        raise ValueError("Value not found!")
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                    if k > 350:
667                        raise NotImplementedError("Maxima is faster here...")
668                    for j in xrange(0,floor(k/2)+1):
669                        help3 = 2*expression.Expression.expand(x*help2) - help1
670                        help1 = help2
671                        help2 = help3
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                if k > 150:
743                    raise NotImplementedError("Maxima is faster here!")
744                for j in xrange(0,k-1):
745                    help3 = 2*expression.Expression.expand(x*help2) - help1
746                    help1 = help2
747                    help2 = help3
748            else:
749                for j in xrange(0,k-1):
750                    help3 = 2*x*help2 - help1
751                    help1 = help2
752                    help2 = help3
753                   
754            return help3   
755   
756    def _maxima_init_evaled_(self, *args):
757        """
758        Uses 
759        EXAMPLES::
760        sage: x = PolynomialRing(QQ, 'x').gen()
761        sage: chebyshev_T(2,x)
762        2*x^2 - 1
763        """
764        n = args[0]
765        x = args[1]
766        return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
767
768       
769    def _evalf_(self, *args,**kwds):
770        """
771        Evals chebyshev_U
772        numerically with mpmath.
773        EXAMPLES::
774            sage: chebyshev_U(10,3).n(75)
775            4.661117900000000000000e7
776        """
777        try:
778            step_parent = kwds['parent']
779        except KeyError:
780            step_parent = parent(args[-1])
781
782        try:
783            precision = step_parent.prec()
784        except AttributeError:
785            precision = mpmath.mp.prec
786
787        return mpmath.call(mpmath.chebyu,args[0],args[-1],prec = precision)
788
789    def _eval_special_values_(self,*args):
790        """
791        Values known at the boundary.
792        EXAMPLES:
793       
794            sage: var('n')
795            n
796            sage: chebyshev_U(n,1)
797            n + 1
798            sage: chebyshev_U(n,-1)
799            (n + 1)*(-1)^n
800        """
801        if args[-1] == 1:
802            return (args[0]+1)
803       
804        if args[-1] == -1:
805            return (-1)**args[0]*(args[0]+1)
806
807        if (args[-1] == 0):
808            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
809
810        raise ValueError("Value not found")
811
812    def _derivative_(self, *args, **kwds):
813        """
814        Returns the derivative of chebyshev_U in form of the chebyshev
815        Polynomials of the first and second kind
816       
817        EXAMPLES::
818            sage: var('k')
819            k
820            sage: derivative(chebyshev_U(k,x),x)
821            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
822            sage: derivative(chebyshev_U(3,x),x)
823            24*x^2 - 4
824            sage: derivative(chebyshev_U(k,x),k)
825            Traceback (most recent call last):
826            ...
827            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
828 
829        """
830        diff_param = kwds['diff_param']
831        if diff_param == 0: 
832            raise NotImplementedError(
833"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
834        else:
835            return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
836                    chebyshev_U(args[0],args[1]))/(args[1]**2-1)
837       
838chebyshev_U = Func_chebyshev_U()
839
840class Func_gegenbauer(OrthogonalPolynomial):
841    """
842    Returns the ultraspherical (or Gegenbauer) polynomial.
843   
844    REFERENCE:
845
846    - AS 22.5.27
847   
848    EXAMPLES::
849   
850        sage: x = PolynomialRing(QQ, 'x').gen()
851        sage: ultraspherical(2,3/2,x)
852        15/2*x^2 - 3/2
853        sage: ultraspherical(2,1/2,x)
854        3/2*x^2 - 1/2
855        sage: ultraspherical(1,1,x)
856        2*x     
857        sage: t = PolynomialRing(RationalField(),"t").gen()
858        sage: gegenbauer(3,2,t)
859        32*t^3 - 12*t
860    """   
861    def __init__(self):
862        OrthogonalPolynomial.__init__(self,"gegenbauer",nargs = 3,conversions =dict(maxima='ultraspherical',mathematica='GegenbauerC'))
863
864    def _clenshaw_method_(self,*args):
865        """
866        Clenshaw method for gegenbauer poly (means use the recursion...)
867        This is much faster for numerical evaluation than maxima!
868        See A.S. 227 (p. 782) for details for the recurions
869        """
870        k = args[0]
871        x = args[-1]
872        alpha = args[1]
873
874        if is_Expression(alpha) or abs(alpha) > numpy.finfo(float).eps:
875            alpha_zero = False
876        else:
877            alpha_zero = True
878
879        if k == 0:
880            return 1
881        elif k == 1:
882            if not alpha_zero:
883                return 2*x*alpha
884            else:
885                return 2*x # It seems that maxima evals this the wrong way! (see A.S. 22.4 (p.777))
886        else:
887            help1 = 1
888            if alpha_zero:
889                help2 = 2*x
890            else:
891                help2 = 2*x*alpha
892
893            help3 = 0
894            if is_Expression(x):
895                for j in xrange(0,floor(k/2)+1):
896                    help3 = help3 + (-1)**j*gamma(alpha+k-j)/factorial(j)/factorial(k-2*j)*(2*x)**(k-2*j)
897            else:
898                for j in xrange(1,k):
899                    help3 = 2*(j+alpha)*x*help2 - (j+2*alpha-1)*help1
900                    help3 = help3/(j+1)
901                    help1 = help2
902                    help2 = help3
903                   
904            return help3   
905
906    def _maxima_init_evaled_(self, *args):
907        """
908        Uses 
909        EXAMPLES::
910        sage: x = PolynomialRing(QQ, 'x').gen()
911        sage: ultraspherical(2,1/2,x)
912        3/2*x^2 - 1/2
913        """
914        n = args[0]
915        a = args[1]
916        x = args[2]
917        return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
918
919
920    def _evalf_(self, *args,**kwds):
921        """
922        Evals gegenbauer
923        numerically with mpmath.
924        EXAMPLES::
925            sage: gegenbauer(10,2,3.).n(54)
926            5.25360702000000e8
927        """
928        try:
929            step_parent = kwds['parent']
930        except KeyError:
931            step_parent = parent(args[-1])
932
933        try:
934            precision = step_parent.prec()
935        except AttributeError:
936            precision = mpmath.mp.prec
937
938        return mpmath.call(mpmath.gegenbauer,args[0],args[1],args[-1],prec = precision)
939
940    def _eval_special_values_(self,*args):
941        """
942        Special values known. A.S. 22.4 (p.777)
943        EXAMPLES:
944       
945            sage: var('n')
946            n
947            sage: chebyshev_U(n,1)
948            n + 1
949            sage: chebyshev_U(n,-1)
950            (n + 1)*(-1)^n
951        """
952        if args[1] == 0 and args[0] != 0:
953            return args[0]*chebyshev_T(args[0],args[-1])/2
954
955        if args[1] == 1/2:
956            return legendre_P(args[0],args[-1])
957
958        if args[1] == 1:
959            return chebyshev_U(args[0],args[-1])
960
961        if args[-1] == 1:
962            return binomial(args[0] + 2*args[1] - 1,args[0])
963       
964        if args[-1] == -1:
965            return (-1)**args[0]*binomial(args[0] + 2*args[1] - 1,args[0])
966
967        if (args[-1] == 0):
968            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2*gamma(args[1]+args[0]/2)/gamma(args[1])/gamma(args[0]/2)
969
970        raise ValueError("Value not found")
971
972gegenbauer = Func_gegenbauer()
973
974class Func_jacobi_P(OrthogonalPolynomial):
975    r"""
976    Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for
977    integers `n > -1` and a and b symbolic or `a > -1`
978    and `b > -1`. The Jacobi polynomials are actually defined
979    for all a and b. However, the Jacobi polynomial weight
980    `(1-x)^a(1+x)^b` isn't integrable for `a \leq -1`
981    or `b \leq -1`.
982   
983    REFERENCE:
984
985    - table on page 789 in AS.
986   
987    EXAMPLES::
988   
989        sage: x = PolynomialRing(QQ, 'x').gen()
990        sage: jacobi_P(2,0,0,x)
991        3/2*x^2 - 1/2
992        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
993        5.009999999999998
994    """
995
996    def __init__(self):
997        OrthogonalPolynomial.__init__(self,"jacobi_P",nargs = 4,
998conversions =dict(maxima='jacobi_p',mathematica='JacobiP'))
999   
1000    def _clenshaw_method_(self,*args):
1001        """
1002        Clenshaw method for jacobi_P (means use the recursion,
1003        or sum formula)
1004        This is much faster for numerical evaluation than maxima!
1005        See A.S. 227 (p. 782) for details for the recurions.
1006        Warning: The clanshaw method for the Jacobi Polynomials
1007        should only used for exact data types, when high orders are
1008        used, due to weak instabilities of the recursion!
1009        """ 
1010        k = args[0]
1011        x = args[-1]
1012        alpha = args[1]
1013        beta = args[2]
1014
1015        if k == 0:
1016            return 1
1017        elif k == 1:
1018            return (alpha-beta + (alpha+beta+2)*x)/2
1019        else:
1020           
1021            if is_Expression(x) or is_Expression(alpha) or is_Expression(beta):
1022                #Here we use the sum formula of jacobi_P it seems this is rather
1023                #optimal for use.
1024                help1 = gamma(alpha+k+1)/factorial(k)/gamma(alpha+beta+k+1)
1025                help2 = 0
1026                for j in xrange(0,k+1):
1027                    help2 = help2 + binomial(k,j)*gamma(alpha+beta+k+j+1)/gamma(alpha+j+1)*((x-1)/2)**j
1028                return help1*help2
1029            else:
1030                help1 = 1
1031                help2 = (alpha-beta + (alpha+beta+2)*x)/2
1032
1033                for j in xrange(1,k):
1034                    a1 = 2*(j+1)*(j+alpha+beta+1)*(2*j+alpha+beta)
1035                    a2 = (2*j+alpha+beta+1)*(alpha**2-beta**2)
1036                    a3 = (2*j+alpha+beta)*(2*j+alpha+beta+1)*(2*j+alpha+beta+2)
1037                    a4 = 2*(j+alpha)*(j+beta)*(2*j+alpha+beta+2)
1038                    help3 = (a2+a3*x)*help2 - a4*help1
1039                    help3 = help3/a1
1040                    help1 = help2
1041                    help2 = help3
1042               
1043            return help3   
1044
1045    def _maxima_init_evaled_(self, *args):
1046        """
1047        Old maxima method.
1048        """
1049        n = args[0]
1050        a = args[1]
1051        b = args[2]
1052        x = args[3]
1053        return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
1054
1055    def _evalf_(self, *args,**kwds):
1056        """
1057        Evals jacobi_P
1058        numerically with mpmath.
1059        EXAMPLES::
1060            sage: jacobi_P(10,2,3,3).n(75)
1061            1.322776620000000000000e8
1062        """
1063        try:
1064            step_parent = kwds['parent']
1065        except KeyError:
1066            step_parent = parent(args[-1])
1067
1068        try:
1069            precision = step_parent.prec()
1070        except AttributeError:
1071            precision = mpmath.mp.prec
1072
1073        return mpmath.call(mpmath.jacobi,args[0],args[1],args[2],args[-1],prec = precision)
1074
1075    def _eval_special_values_(self,*args):
1076        """
1077        Special values known.
1078        EXAMPLES:
1079           
1080            sage: var('n')
1081            n
1082            sage: legendre_P(n,1)
1083            1
1084            sage: legendre_P(n,-1)
1085            (-1)^n
1086        """
1087        if args[-1] == 1:
1088            return binomial(args[0]+args[1],args[0])
1089       
1090        if args[-1] == -1:
1091            return (-1)**args[0]*binomial(args[0]+args[1],args[0])
1092       
1093        if args[1] == 0 and args[2] == 0:
1094            return legendre_P(args[0],args[-1])
1095
1096        if args[1] == -0.5 and args[2] == -0.5:
1097            try:
1098                return binomial(2*args[0],args[0])*chebyshev_T(args[0],args[-1])/4**args[0]
1099            except TypeError:
1100                pass
1101
1102        raise ValueError("Value not found")
1103
1104    def _derivative_(self, *args, **kwds):
1105        """
1106        Returns the derivative of jacobi_P in form of jacobi_polynomials
1107       
1108        EXAMPLES::
1109            sage: var('k a b')
1110            (k, a, b)
1111            sage: derivative(jacobi_P(k,a,b,x),x)
1112            -(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))
1113            sage: derivative(jacobi_P(2,1,3,x),x)
1114            14*x - 7/2
1115            sage: derivative(jacobi_P(k,a,b,x),a)
1116            Traceback (most recent call last):
1117            ...
1118            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
1119 
1120        """
1121        diff_param = kwds['diff_param']
1122        if diff_param in [0,1,2]: 
1123            raise NotImplementedError(
1124"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1125        else:
1126            return (args[0]*(args[1]-args[2]-(2*args[0]+args[1]+args[2])*args[-1])*
1127                    jacobi_P(args[0],args[1],args[2],args[-1])+ 2*(args[0]+args[1])*
1128                    (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)
1129       
1130       
1131jacobi_P = Func_jacobi_P()       
1132
1133class Func_legendre_P(OrthogonalPolynomial):
1134    """
1135    Returns the Legendre polynomial of the first kind..
1136   
1137    REFERENCE:
1138
1139    - AS 22.5.35 page 779.
1140   
1141    EXAMPLES::
1142   
1143        sage: P.<t> = QQ[]
1144        sage: legendre_P(2,t)
1145        3/2*t^2 - 1/2
1146        sage: legendre_P(3, 1.1)
1147        1.67750000000000
1148        sage: legendre_P(3, 1 + I)
1149        7/2*I - 13/2
1150        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
1151        Traceback (most recent call last):
1152        ...
1153        TypeError: cannot coerce arguments: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Integer Ring to Symbolic Ring
1154        sage: legendre_P(3, GF(11)(5))
1155        8
1156    """
1157    def __init__(self):
1158        OrthogonalPolynomial.__init__(self,"legendre_P",nargs = 2,
1159conversions =dict(maxima='legendre_p',mathematica='LegendreP'))
1160
1161    def _clenshaw_method_(self,*args):
1162        """
1163        Clenshaw method for legendre_P (means use the recursion...)
1164        This is much faster for numerical evaluation than maxima!
1165        See A.S. 227 (p. 782) for details for the recurions.
1166        Warning: The clanshaw method for the Legendre Polynomials
1167        should only used for exact data types, when high orders are
1168        used, due to weak instabilities of the recursion!
1169        """
1170        k = args[0]
1171        x = args[-1]
1172       
1173        if k == 0:
1174            return 1
1175        elif k == 1:
1176            return x
1177        else:
1178            help1 = 1
1179            help2 = x
1180            if is_Expression(x):
1181                if k > 100:
1182                    raise NotImplementedError("Maxima is faster here...")
1183               
1184                for j in xrange(1,k):
1185                    help3 = (2*j+1)*expression.Expression.expand(x*help2)- j*help1
1186                    help3 = help3/(j+1)
1187                    help1 = help2
1188                    help2 = help3
1189            else:
1190                for j in xrange(1,k):
1191                    help3 = (2*j+1)*x*help2 - j*help1
1192                    help3 = help3/(j+1)
1193                    help1 = help2
1194                    help2 = help3
1195               
1196            return help3   
1197
1198    def _maxima_init_evaled_(self, *args):
1199        n = args[0]
1200        x = args[1]
1201        return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)), locals={'x':x})
1202
1203    def _evalf_(self, *args,**kwds):
1204        """
1205        Evals legendre_P
1206        numerically with mpmath.
1207        EXAMPLES::
1208            sage: legendre_P(10,3).n(75)
1209            8.097453000000000000000e6
1210        """
1211        try:
1212            step_parent = kwds['parent']
1213        except KeyError:
1214            step_parent = parent(args[-1])
1215
1216        try:
1217            precision = step_parent.prec()
1218        except AttributeError:
1219            precision = mpmath.mp.prec
1220
1221        return mpmath.call(mpmath.legendre,args[0],args[-1],prec = precision)
1222   
1223
1224    def _eval_special_values_(self,*args):
1225        """
1226        Special values known.
1227        EXAMPLES:
1228           
1229            sage: var('n')
1230            n
1231            sage: legendre_P(n,1)
1232            1
1233            sage: legendre_P(n,-1)
1234            (-1)^n
1235        """
1236        if args[-1] == 1:
1237            return 1
1238       
1239        if args[-1] == -1:
1240            return (-1)**args[0]
1241       
1242        if (args[-1] == 0):
1243            try:
1244                return (1+(-1)**args[0])/2*binomial(args[0],args[0]/2)/4**(args[0]/2)
1245            except TypeError:
1246                pass
1247
1248        raise ValueError("Value not found")
1249
1250    def _derivative_(self,*args,**kwds):
1251        """return the derivative of legendre_P in
1252        form of the Legendre Polynomial.
1253        EXAMPLES::
1254        derivative(legendre_P(n,x),x)
1255        -(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x - 1)^2
1256        derivative(legendre_P(3,x),x)
1257        15/2*x^2 - 3/2
1258        """
1259        diff_param = kwds['diff_param']
1260        if diff_param == 0: 
1261            raise NotImplementedError(
1262"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1263        else:
1264            return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*legendre_P(args[0],args[-1]))/(1-args[-1])**2
1265       
1266
1267legendre_P = Func_legendre_P()     
1268
1269class Func_legendre_Q(OrthogonalPolynomial):
1270    """
1271    Return the chebyshev function of the second kind.
1272   
1273    REFERENCE:
1274   
1275    - A.S. 8 (p. 332)
1276   
1277    EXAMPLES::
1278   
1279        sage: P.<t> = QQ[]
1280        sage: legendre_Q(2, t)
1281        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
1282        sage: legendre_Q(3, 0.5)
1283        -0.198654771479482
1284        sage: legendre_Q(4, 2)
1285        443/16*I*pi + 443/16*log(3) - 365/12
1286        sage: legendre_Q(4, 2.0)
1287        0.00116107583162324 + 86.9828465962674*I
1288
1289    """
1290    def __init__(self):
1291        OrthogonalPolynomial.__init__(self,"legendre_Q",nargs = 2,
1292conversions =dict(maxima='legendre_q',mathematica='LegendreQ'))
1293
1294    def _maxima_init_evaled_(self, *args):
1295        """
1296        Maxima seems just fine for legendre Q. So we use it here!
1297        """
1298        n = args[0]
1299        x = args[1]
1300        return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
1301
1302    def _clenshaw_method_(self,*args):
1303        """
1304        Clenshaw method for legendre_q (means use the recursion...)
1305        This is much faster for numerical evaluation than maxima!
1306        See A.S. 8.5.3 (p. 334) for details for the recurions.
1307        Warning: The clanshaw method for the Legendre fUNCTIONS
1308        should only used for exact data types, when high orders are
1309        used, due to weak instabilities of the recursion!
1310        """
1311        k = args[0]
1312        x = args[-1]
1313       
1314        if k == 0:
1315            return ln((1+x)/(1-x))/2
1316        elif k == 1:
1317            return x/2*ln((1+x)/(1-x))-1
1318        else:
1319            if is_Expression(x):
1320                raise NotImplementedError("Maxima works fine here!") 
1321                #it seems that the old method just works fine here...
1322                #raise NotImplementedError("clenshaw does not work well...")
1323            else:
1324                help1 = ln((1+x)/(1-x))/2
1325                help2 = x/2*ln((1+x)/(1-x))-1
1326
1327                for j in xrange(1,k):
1328                    help3 = (2*j+1)*x*help2 - j*help1
1329                    help3 = help3/(j+1)
1330                    help1 = help2
1331                    help2 = help3
1332               
1333            return help3   
1334
1335    def _eval_special_values_(self,*args):
1336        """
1337        Special values known.
1338        EXAMPLES:
1339           
1340            sage: var('n')
1341            n
1342            sage: legendre_P(n,1)
1343            1
1344            sage: legendre_P(n,-1)
1345            (-1)^n
1346        """
1347        if args[-1] == 1:
1348            return NaN
1349       
1350        if args[-1] == -1:
1351            return NaN
1352       
1353        if (args[-1] == 0):
1354            if is_Expression(args[0]): 
1355                try:
1356                    return -(SR.pi()**(QQ(1)/QQ(2)))/2*sin(SR.pi()/2*args[0])*gamma((args[0]+1)/2)*gamma(args[0]/2 + 1)
1357                except TypeError:
1358                    pass
1359            else:
1360                return -(math.pi**(QQ(1)/QQ(2)))/2*sin(math.pi/2*args[0])*gamma((args[0]+1)/2)*gamma(args[0]/2. + 1)
1361
1362        raise ValueError("Value not found")
1363
1364    def _evalf_(self, *args,**kwds):
1365        """
1366        Evals legendre_Q
1367        numerically with mpmath.
1368        Function does not work yet...
1369        EXAMPLES::
1370            sage: legendre_P(10,3).n(75)
1371            8.097453000000000000000e6
1372        """
1373        raise AttributeError("mpmath function does not work correctly!")
1374
1375    def _derivative_(self,*args,**kwds):
1376        """return the derivative of legendre_P in
1377        form of the Legendre Polynomial.
1378        EXAMPLES::
1379        derivative(legendre_P(n,x),x)
1380        -(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x - 1)^2
1381        derivative(legendre_P(3,x),x)
1382        15/2*x^2 - 3/2
1383        """
1384        diff_param = kwds['diff_param']
1385        if diff_param == 0: 
1386            raise NotImplementedError(
1387"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1388        else:
1389            raise NotImplementedError("Not implemented yet")
1390            return (gen_legendre_Q(args[0],1,args[-1]))/sqrt(1-args[-1]**2)
1391
1392legendre_Q = Func_legendre_Q()
1393
1394def legendre_Q_old(n,x):
1395    """
1396    Returns the Legendre function of the second kind for integers
1397    `n>-1`. (old version)
1398   
1399    Computed using Maxima.
1400   
1401    EXAMPLES::
1402   
1403        sage: P.<t> = QQ[]
1404        sage: legendre_Q(2, t)
1405        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
1406        sage: legendre_Q(3, 0.5)
1407        -0.198654771479482
1408        sage: legendre_Q(4, 2)
1409        443/16*I*pi + 443/16*log(3) - 365/12
1410        sage: legendre_Q(4, 2.0)
1411        0.00116107583162324 + 86.9828465962674*I
1412    """
1413    _init()
1414    return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
1415
1416
1417
1418def gen_laguerre(n,a,x):
1419    """
1420    Returns the generalized Laguerre polynomial for integers `n > -1`.
1421    Typically, a = 1/2 or a = -1/2.
1422   
1423    REFERENCE:
1424
1425    - table on page 789 in AS.
1426   
1427    EXAMPLES::
1428   
1429        sage: x = PolynomialRing(QQ, 'x').gen()
1430        sage: gen_laguerre(2,1,x)
1431        1/2*x^2 - 3*x + 3
1432        sage: gen_laguerre(2,1/2,x)
1433        1/2*x^2 - 5/2*x + 15/8
1434        sage: gen_laguerre(2,-1/2,x)
1435        1/2*x^2 - 3/2*x + 3/8
1436        sage: gen_laguerre(2,0,x)
1437        1/2*x^2 - 2*x + 1
1438        sage: gen_laguerre(3,0,x)
1439        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1440    """
1441    from sage.functions.all import sqrt
1442    _init()
1443    return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
1444
1445def gen_legendre_P(n,m,x):
1446    r"""
1447    Returns the generalized (or associated) Legendre function of the
1448    first kind for integers `n > -1, m > -1`.
1449   
1450    The awkward code for when m is odd and 1 results from the fact that
1451    Maxima is happy with, for example, `(1 - t^2)^3/2`, but
1452    Sage is not. For these cases the function is computed from the
1453    (m-1)-case using one of the recursions satisfied by the Legendre
1454    functions.
1455   
1456    REFERENCE:
1457
1458    - Gradshteyn and Ryzhik 8.706 page 1000.
1459   
1460    EXAMPLES::
1461   
1462        sage: P.<t> = QQ[]
1463        sage: gen_legendre_P(2, 0, t)
1464        3/2*t^2 - 1/2
1465        sage: gen_legendre_P(2, 0, t) - legendre_P(2, t)
1466        0
1467        sage: gen_legendre_P(3, 1, t)
1468        -3/2*sqrt(-t^2 + 1)*(5*t^2 - 1)
1469        sage: gen_legendre_P(4, 3, t)
1470        105*sqrt(-t^2 + 1)*(t^2 - 1)*t
1471        sage: gen_legendre_P(7, 3, I).expand()
1472        -16695*sqrt(2)
1473        sage: gen_legendre_P(4, 1, 2.5)
1474        -583.562373654533*I
1475    """
1476    from sage.functions.all import sqrt
1477    _init()
1478    if m.mod(2).is_zero() or m.is_one():
1479        return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
1480    else:
1481        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))
1482
1483def gen_legendre_Q(n,m,x):
1484    """
1485    Returns the generalized (or associated) Legendre function of the
1486    second kind for integers `n>-1`, `m>-1`.
1487   
1488    Maxima restricts m = n. Hence the cases m n are computed using the
1489    same recursion used for gen_legendre_P(n,m,x) when m is odd and
1490    1.
1491   
1492    EXAMPLES::
1493   
1494        sage: P.<t> = QQ[]
1495        sage: gen_legendre_Q(2,0,t)
1496        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
1497        sage: gen_legendre_Q(2,0,t) - legendre_Q(2, t)
1498        0
1499        sage: gen_legendre_Q(3,1,0.5)
1500        2.49185259170895
1501        sage: gen_legendre_Q(0, 1, x)
1502        -1/sqrt(-x^2 + 1)
1503        sage: gen_legendre_Q(2, 4, x).factor()
1504        48*x/((x - 1)^2*(x + 1)^2)
1505    """
1506    from sage.functions.all import sqrt
1507    if m <= n:
1508        _init()
1509        return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
1510    if m == n + 1 or n == 0:
1511        if m.mod(2).is_zero():
1512            denom = (1 - x**2)**(m/2)
1513        else:
1514            denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)   
1515        if m == n + 1:
1516            return (-1)**m*(m-1).factorial()*2**n/denom
1517        else:
1518            return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
1519    else:
1520        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)
1521
1522def hermite(n,x):
1523    """
1524    Returns the Hermite polynomial for integers `n > -1`.
1525   
1526    REFERENCE:
1527
1528    - AS 22.5.40 and 22.5.41, page 779.
1529   
1530    EXAMPLES::
1531   
1532        sage: x = PolynomialRing(QQ, 'x').gen()
1533        sage: hermite(2,x)
1534        4*x^2 - 2
1535        sage: hermite(3,x)
1536        8*x^3 - 12*x
1537        sage: hermite(3,2)
1538        40
1539        sage: S.<y> = PolynomialRing(RR)
1540        sage: hermite(3,y)
1541        8.00000000000000*y^3 - 12.0000000000000*y
1542        sage: R.<x,y> = QQ[]
1543        sage: hermite(3,y^2)
1544        8*y^6 - 12*y^2
1545        sage: w = var('w')
1546        sage: hermite(3,2*w)
1547        8*(8*w^2 - 3)*w
1548    """
1549    _init()
1550    return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
1551
1552def jacobi_P_old(n,a,b,x):
1553    r"""
1554    Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for
1555    integers `n > -1` and a and b symbolic or `a > -1`
1556    and `b > -1`. The Jacobi polynomials are actually defined
1557    for all a and b. However, the Jacobi polynomial weight
1558    `(1-x)^a(1+x)^b` isn't integrable for `a \leq -1`
1559    or `b \leq -1`.
1560   
1561    REFERENCE:
1562
1563    - table on page 789 in AS.
1564   
1565    EXAMPLES::
1566   
1567        sage: x = PolynomialRing(QQ, 'x').gen()
1568        sage: jacobi_P(2,0,0,x)
1569        3/2*x^2 - 1/2
1570        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
1571        5.009999999999998
1572    """
1573    _init()
1574    return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
1575
1576def laguerre(n,x):
1577    """
1578    Returns the Laguerre polynomial for integers `n > -1`.
1579   
1580    REFERENCE:
1581
1582    - AS 22.5.16, page 778 and AS page 789.
1583   
1584    EXAMPLES::
1585   
1586        sage: x = PolynomialRing(QQ, 'x').gen()
1587        sage: laguerre(2,x)
1588        1/2*x^2 - 2*x + 1
1589        sage: laguerre(3,x)
1590        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1591        sage: laguerre(2,2)
1592        -1
1593    """
1594    _init()
1595    return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
1596
1597def ultraspherical(n,a,x):
1598    """
1599    Returns the ultraspherical (or Gegenbauer) polynomial for integers
1600    `n > -1`.
1601   
1602    Computed using Maxima.
1603   
1604    REFERENCE:
1605
1606    - AS 22.5.27
1607   
1608    EXAMPLES::
1609   
1610        sage: x = PolynomialRing(QQ, 'x').gen()
1611        sage: ultraspherical(2,3/2,x)
1612        15/2*x^2 - 3/2
1613        sage: ultraspherical(2,1/2,x)
1614        3/2*x^2 - 1/2
1615        sage: ultraspherical(1,1,x)
1616        2*x     
1617        sage: t = PolynomialRing(RationalField(),"t").gen()
1618        sage: gegenbauer(3,2,t)
1619        32*t^3 - 12*t
1620    """
1621    _init()
1622    return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
1623
1624