Ticket #9706: orthogonal_polys.2.py

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

Newer version, with legendre_P, and faster evaluation of symbolic expressions

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
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
330import sage.symbolic.expression as expression
331
332_done = False
333def _init():
334    """
335    Internal function which checks if Maxima has loaded the
336    "orthopoly" package.  All functions using this in this
337    file should call this function first.
338   
339    TEST:
340
341    The global starts ``False``::
342
343        sage: sage.functions.orthogonal_polys._done
344        False
345
346    Then after using one of these functions, it changes:: (The value is now
347    False for chebyshev_T because chebyshev_T uses clenshaw method instead...)
348
349        sage: from sage.functions.orthogonal_polys import chebyshev_T
350        sage: chebyshev_T(2,x)
351        2*x^2 - 1
352        sage: sage.functions.orthogonal_polys._done
353        False
354        sage: gegenbauer(1,2,3)
355        12
356        sage: sage.functions.orthogonal_polys._done
357        True
358   
359    Note that because here we use a Pynac variable ``x``,
360    the representation of the function is different from
361    its actual doctest, where a polynomial indeterminate
362    ``x`` is used.
363    """
364    global _done
365    if _done:
366        return
367    maxima.eval('load("orthopoly");')
368    # TODO -- make it possible to use the intervals returned
369    # instead of just discarding this info!
370    maxima.eval('orthopoly_returns_intervals:false;')
371    _done = True
372
373#load /home/maldun/sage/sage-4.5.1/devel/sage-ortho/sage/functions/orthogonal_polys.py
374
375class OrthogonalPolynomial(BuiltinFunction):
376    """
377    Base Class for Orthogonal Polynomials. The evaluation as a polynomial
378    is done via maxima due performance reasons. Therefore the internal name
379    in maxima maxima_name has to be declared.
380    Convention: The first argument is always the order of the polynomial,
381    he last one is always the value x where the polynomial is evaluated.
382
383    """
384    def __init__(self, name, maxima_name, nargs = 2, latex_name = None, conversions = None):
385        self._maxima_name = maxima_name
386        BuiltinFunction.__init__(self, name = name, nargs = nargs, latex_name = latex_name, conversions = conversions)
387
388    def _maxima_init_evaled_(self, *args):
389        """
390        Returns a string which represents this function evaluated at
391        *args* in Maxima.
392        Remark: This function is stolen from the class MaximaFunction,
393        from sage.functions.special.py
394        Comment: The usefulness is in question....
395
396        EXAMPLES::
397
398            sage: chebyshev_T(3,x)
399            4*x^3 - 3*x
400        """
401        args_maxima = []
402        for a in args:
403            if isinstance(a, str):
404                args_maxima.append(a)
405            elif hasattr(a, '_maxima_init_'):
406                args_maxima.append(a._maxima_init_())
407            else:
408                args_maxima.append(str(a))
409        return "%s(%s)"%(self._maxima_name, ', '.join(args_maxima))
410
411    def _clenshaw_method_(self,*args):
412        """
413        The Clenshaw method uses the three term recursion of the polynomial,
414        or explicit formulas instead of maxima to evaluate the polynomial
415        efficiently, if the x argument is not a symbolic expression.
416        The name comes from the Clenshaw algorithm for fast evaluation of
417        polynomialsin chebyshev form.
418       
419        comparison: Maxima and Clenshaw algorithm for non-symbolic evaluation:
420         #sage: time chebyshev_T(50,10) #clenshaw
421         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
422         #Wall time: 0.00 s
423         #49656746733678312490954442369580252421769338391329426325400124999
424         #sage: time sage.functions.orthogonal_polys.chebyshev_T(50,10)
425         #maxima
426         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
427         #Wall time: 0.05 s
428         #49656746733678312490954442369580252421769338391329426325400124999
429
430         #sage: time chebyshev_T(500,10); #clenshaw
431         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
432         #Wall time: 0.00 s
433         #sage: time sage.functions.orthogonal_polys.chebyshev_T(500,10);
434         #maxima
435         #CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s
436         #Wall time: 0.77 s
437        """
438        raise NotImplementedError(
439            "No recursive calculation of values implemented (yet)!")
440
441    def _eval_special_values_(self,*args):
442        """
443        Evals the polynomial explicitly for special values.
444        EXAMPLES:
445           
446            sage: var('n')
447            n
448            sage: chebyshev_T(n,-1)
449            (-1)^n
450        """
451        raise ValueError("No special values known!")
452
453    def _eval_(self, *args):
454        """
455       
456        The symbolic evaluation is done with maxima, because the evaluation of
457        the Polynomial representation seems to be quite clever.
458        For the fast numerical evaluation an other method should be used...
459        Therefore I suggest Clenshaw's algorithm, which uses the rekursion!
460        The function also checks for special values, and if
461        the order is an integer and in range!
462       
463        Update: I take back my argument about maxima...
464        it seems that just adding the expand command to
465        the for loop in the Clanshaw algorithm speeds up
466        the symbolic evaluation enormously.
467        Now the maxima option stays in the class, because
468        it could be good for more complicated polynomials.
469        But perhaps I will throw it out later...
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            0
500            sage: chebyshev_T(7.5,x)
501            Traceback (most recent call last):
502            ...
503            TypeError: Order is not an Integer!
504
505       
506           
507           
508        """
509        if not is_Expression(args[0]):
510            if args[0] != floor(args[0]):
511                raise TypeError("Order is not an Integer!")
512       
513            if args[0] < 0:
514                return 0 #raise TypeError("Order is not valid!")
515
516        try:
517            return self._eval_special_values_(*args)
518        except ValueError:
519            pass
520
521         
522        if not is_Expression(args[0]):
523            try:
524                return self._clenshaw_method_(*args)
525            except NotImplementedError:
526                pass
527
528        _init()
529        try:
530            s = maxima(self._maxima_init_evaled_(*args))
531        except TypeError:
532            return None
533        if self._maxima_name in repr(s):
534            return None
535        else:
536            return s.sage()
537   
538           
539        #TODO: numpy_eval with help of the new scipy version!!!!
540        #Reason scipy 0.8 supports stable and fast numerical evaluation
541        #of ortho polys
542
543       
544class Func_chebyshev_T(OrthogonalPolynomial): 
545
546    """
547    Class for the Chebyshev polynomial of the first kind for integers
548    `n>-1`.
549   
550    REFERENCE:
551
552    - AS 22.5.31 page 778 and AS 6.1.22 page 256.
553   
554    EXAMPLES::
555   
556       sage: chebyshev_T(3,x)
557       4*x^3 - 3*x
558       sage: var('k')
559       k
560       sage: test = chebyshev_T(k,x)
561       sage: test
562       chebyshev_T(k, x)
563    """
564   
565    def __init__(self):
566        OrthogonalPolynomial.__init__(self,"chebyshev_T","chebyshev_t",nargs = 2,
567conversions =dict(maxima='chebyshev_t',mathematica='ChebyshevT'))
568   
569    def _eval_special_values_(self,*args):
570        """
571        Values known for special values of x.
572        For details see A.S. 22.4 (p. 777)
573       
574        EXAMPLES:
575           
576            sage: var('n')
577            n
578            sage: chebyshev_T(n,1)
579            1
580            sage: chebyshev_T(n,-1)
581            (-1)^n
582        """
583        if args[-1] == 1:
584            return 1
585       
586        if args[-1] == -1:
587            return (-1)**args[0]
588
589        if (args[-1] == 0):
590            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
591
592        raise ValueError("Value not found!")
593   
594    def _clenshaw_method_(self,*args):
595            """
596            Clenshaw method for chebyshev_T (means use recursions in this case)
597            This is much faster for numerical evaluation than maxima!
598            See A.S. 227 (p. 782) for details for the recurions
599            """
600
601            k = args[0]
602            x = args[1]
603
604            if k == 0:
605                return 1
606            elif k == 1:
607                return x
608            else:
609                #TODO: When evaluation of Symbolic Expressions works better
610                #use these explicit formulas instead!
611                #if -1 <= x <= 1:
612                #    return cos(k*acos(x))
613                #elif 1 < x:
614                #    return cosh(k*acosh(x))
615                #else: # x < -1
616                #    return (-1)**(k%2)*cosh(k*acosh(-x))
617               
618                help1 = 1
619                help2 = x
620                if is_Expression(x):
621                    for j in xrange(0,k-1):
622                        help3 = expression.Expression.expand(2*x*help2 - help1)
623                        help1 = help2
624                        help2 = help3
625                else:
626                     for j in xrange(0,k-1):
627                        help3 = 2*x*help2 - help1
628                        help1 = help2
629                        help2 = help3
630               
631                return help3
632
633    def _evalf_(self, *args,**kwds):
634        """
635        Evals chebyshev_T numerically with mpmath.
636        EXAMPLES::
637            sage: chebyshev_T(10,3).n(75)
638            2.261953700000000000000e7
639        """
640        parent = kwds['parent']
641        return mpmath.call(mpmath.chebyt,args[0],args[-1],prec = parent.prec())
642
643    def _derivative_(self, *args, **kwds):
644        """
645        Returns the derivative of chebyshev_T in form of the chebyshev Polynomial
646        of the second kind chebyshev_U
647        EXAMPLES::
648            sage: var('k')
649            k
650            sage: derivative(chebyshev_T(k,x),x)
651            k*chebyshev_U(k - 1, x)
652            sage: derivative(chebyshev_T(3,x),x)
653            12*x^2 - 3
654            sage: derivative(chebyshev_T(k,x),k)
655            Traceback (most recent call last):
656            ...
657            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
658 
659        """
660        diff_param = kwds['diff_param']
661        if diff_param == 0: 
662            raise NotImplementedError(
663"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
664        else:
665            return args[0]*chebyshev_U(args[0]-1,args[1])
666
667chebyshev_T = Func_chebyshev_T()
668
669class Func_chebyshev_U(OrthogonalPolynomial): 
670   
671    """
672    Class for the Chebyshev polynomial of the second kind for integers `n>=-1`.
673   
674    REFERENCE:
675
676    - AS, 22.8.3 page 783 and AS 6.1.22 page 256.
677   
678    EXAMPLES::
679   
680        sage: chebyshev_U(3,x)
681        8*x^3 - 4*x
682    """
683
684    def __init__(self):
685        OrthogonalPolynomial.__init__(self,"chebyshev_U","chebyshev_u",
686nargs = 2,conversions =dict(maxima='chebyshev_u',mathematica='ChebyshevU'))
687       
688    def _clenshaw_method_(self,*args):
689            """
690            Clenshaw method for chebyshev_U (means use the recursion...)
691            This is much faster for numerical evaluation than maxima!
692            See A.S. 227 (p. 782) for details for the recurions
693            """
694   
695            k = args[0]
696            x = args[1]
697           
698            if k == 0:
699                return 1
700            elif k == 1:
701                return 2*x
702            else:
703                help1 = 1
704                help2 = 2*x
705                if is_Expression(x):
706                    for j in xrange(0,k-1):
707                        help3 = expression.Expression.expand(2*x*help2 - help1)
708                        help1 = help2
709                        help2 = help3
710                else:
711                    for j in xrange(0,k-1):
712                        help3 = 2*x*help2 - help1
713                        help1 = help2
714                        help2 = help3
715
716                return help3   
717     
718    def _eval_special_values_(self,*args):
719        """
720        Values known at the boundary.
721        EXAMPLES:
722
723            sage: var('n')
724            n
725            sage: chebyshev_U(n,1)
726            n + 1
727            sage: chebyshev_U(n,-1)
728            (n + 1)*(-1)^n
729        """
730        if args[-1] == 1:
731            return (args[0]+1)
732       
733        if args[-1] == -1:
734            return (-1)**args[0]*(args[0]+1)
735
736        if (args[-1] == 0):
737            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
738
739        raise ValueError("Value not found")
740
741    def _evalf_(self, *args,**kwds):
742        """
743        Evals chebyshev_U numerically with mpmath.
744        EXAMPLES::
745            sage: chebyshev_U(10,3).n(75)
746            4.661117900000000000000e7
747        """
748        parent = kwds['parent']
749        return mpmath.call(mpmath.chebyu,args[0],args[-1],prec = parent.prec())
750
751    def _derivative_(self, *args, **kwds):
752        """
753        Returns the derivative of chebyshev_U in form of the chebyshev
754        Polynomials of the first and second kind
755       
756        EXAMPLES::
757            sage: var('k')
758            k
759            sage: derivative(chebyshev_U(k,x),x)
760            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
761            sage: derivative(chebyshev_U(3,x),x)
762            24*x^2 - 4
763            sage: derivative(chebyshev_U(k,x),k)
764            Traceback (most recent call last):
765            ...
766            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
767 
768        """
769        diff_param = kwds['diff_param']
770        if diff_param == 0: 
771            raise NotImplementedError(
772"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
773        else:
774            return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
775                    chebyshev_U(args[0],args[1]))/(args[1]**2-1)
776       
777chebyshev_U = Func_chebyshev_U()
778
779class Func_legendre_P(OrthogonalPolynomial):
780    """
781    Returns the Legendre polynomial of the first kind for integers
782    `n >= -1`.
783   
784    REFERENCE:
785
786    - AS 22.5.35 page 779.
787   
788    EXAMPLES::
789   
790        sage: P.<t> = QQ[]
791        sage: legendre_P(2,t)
792        3/2*t^2 - 1/2
793        sage: legendre_P(3, 1.1)
794        1.67750000000000
795        sage: legendre_P(3, 1 + I)
796        7/2*I - 13/2
797        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
798        Traceback (most recent call last):
799        ...
800        TypeError: cannot coerce arguments: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Integer Ring to Symbolic Ring
801        sage: legendre_P(3, GF(11)(5))
802        8
803    """
804    def __init__(self):
805        OrthogonalPolynomial.__init__(self,"legendre_P","legendre_p",nargs = 2,conversions =dict(maxima='legendre_p',mathematica='LegendreP'))
806
807    def _clenshaw_method_(self,*args):
808        """
809        Clenshaw method for legendre_P (means use the recursion...)
810        This is much faster for numerical evaluation than maxima!
811        See A.S. 227 (p. 782) for details for the recurions.
812        Warning: The clanshaw method for the Legendre Polynomials
813        should only used for exact data types, when high orders are
814        used, due to weak instabilities of the recursion!
815        """
816        k = args[0]
817        x = args[-1]
818       
819        if k == 0:
820            return 1
821        elif k == 1:
822            return x
823        else:
824            help1 = 1
825            help2 = x
826            if is_Expression(x):
827                for j in xrange(1,k):
828                    help3 = (2*j+1)*x*help2 - j*help1
829                    help3 = expression.Expression.expand(help3/(j+1))
830                    help1 = help2
831                    help2 = help3
832            else:
833                for j in xrange(1,k):
834                    help3 = (2*j+1)*x*help2 - j*help1
835                    help3 = help3/(j+1)
836                    help1 = help2
837                    help2 = help3
838               
839            return help3   
840
841    def _eval_special_values_(self,*args):
842        """
843        Special values known.
844        EXAMPLES:
845           
846            sage: var('n')
847            n
848            sage: legendre_P(n,1)
849            1
850            sage: legendre_P(n,-1)
851            (-1)^n
852        """
853        if args[-1] == 1:
854            return 1
855       
856        if args[-1] == -1:
857            return (-1)**args[0]
858       
859        if (args[-1] == 0):
860            try:
861                return (1+(-1)**args[0])/2*binomial(args[0],args[0]/2)/4**(args[0]/2)
862            except TypeError:
863                pass
864
865        raise ValueError("Value not found")
866
867    def _evalf_(self, *args,**kwds):
868        """
869        Evals Legendre_P numerically with mpmath.
870        EXAMPLES::
871            sage: legendre_P(10,5).n(53)
872            1.60047267700000e9
873        """
874        parent = kwds['parent']
875        return mpmath.call(mpmath.legendre,args[0],args[-1],prec = parent.prec())
876
877    def _derivative_(self,*args,**kwds):
878        """return the derivative of legendre_P in
879        form of the Legendre Polynomial.
880        EXAMPLES::
881        derivative(legendre_P(n,x),x)
882        -(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x - 1)^2
883        derivative(legendre_P(3,x),x)
884        15/2*x^2 - 3/2
885        """
886        diff_param = kwds['diff_param']
887        if diff_param == 0: 
888            raise NotImplementedError(
889"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
890        else:
891            return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*legendre_P(args[0],args[-1]))/(1-args[-1])**2
892       
893
894legendre_P = Func_legendre_P()     
895
896def gen_laguerre(n,a,x):
897    """
898    Returns the generalized Laguerre polynomial for integers `n > -1`.
899    Typically, a = 1/2 or a = -1/2.
900   
901    REFERENCE:
902
903    - table on page 789 in AS.
904   
905    EXAMPLES::
906   
907        sage: x = PolynomialRing(QQ, 'x').gen()
908        sage: gen_laguerre(2,1,x)
909        1/2*x^2 - 3*x + 3
910        sage: gen_laguerre(2,1/2,x)
911        1/2*x^2 - 5/2*x + 15/8
912        sage: gen_laguerre(2,-1/2,x)
913        1/2*x^2 - 3/2*x + 3/8
914        sage: gen_laguerre(2,0,x)
915        1/2*x^2 - 2*x + 1
916        sage: gen_laguerre(3,0,x)
917        -1/6*x^3 + 3/2*x^2 - 3*x + 1
918    """
919    from sage.functions.all import sqrt
920    _init()
921    return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
922
923def gen_legendre_P(n,m,x):
924    r"""
925    Returns the generalized (or associated) Legendre function of the
926    first kind for integers `n > -1, m > -1`.
927   
928    The awkward code for when m is odd and 1 results from the fact that
929    Maxima is happy with, for example, `(1 - t^2)^3/2`, but
930    Sage is not. For these cases the function is computed from the
931    (m-1)-case using one of the recursions satisfied by the Legendre
932    functions.
933   
934    REFERENCE:
935
936    - Gradshteyn and Ryzhik 8.706 page 1000.
937   
938    EXAMPLES::
939   
940        sage: P.<t> = QQ[]
941        sage: gen_legendre_P(2, 0, t)
942        3/2*t^2 - 1/2
943        sage: gen_legendre_P(2, 0, t) - legendre_P(2, t)
944        0
945        sage: gen_legendre_P(3, 1, t)
946        -3/2*sqrt(-t^2 + 1)*(5*t^2 - 1)
947        sage: gen_legendre_P(4, 3, t)
948        105*sqrt(-t^2 + 1)*(t^2 - 1)*t
949        sage: gen_legendre_P(7, 3, I).expand()
950        -16695*sqrt(2)
951        sage: gen_legendre_P(4, 1, 2.5)
952        -583.562373654533*I
953    """
954    from sage.functions.all import sqrt
955    _init()
956    if m.mod(2).is_zero() or m.is_one():
957        return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
958    else:
959        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))
960
961def gen_legendre_Q(n,m,x):
962    """
963    Returns the generalized (or associated) Legendre function of the
964    second kind for integers `n>-1`, `m>-1`.
965   
966    Maxima restricts m = n. Hence the cases m n are computed using the
967    same recursion used for gen_legendre_P(n,m,x) when m is odd and
968    1.
969   
970    EXAMPLES::
971   
972        sage: P.<t> = QQ[]
973        sage: gen_legendre_Q(2,0,t)
974        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
975        sage: gen_legendre_Q(2,0,t) - legendre_Q(2, t)
976        0
977        sage: gen_legendre_Q(3,1,0.5)
978        2.49185259170895
979        sage: gen_legendre_Q(0, 1, x)
980        -1/sqrt(-x^2 + 1)
981        sage: gen_legendre_Q(2, 4, x).factor()
982        48*x/((x - 1)^2*(x + 1)^2)
983    """
984    from sage.functions.all import sqrt
985    if m <= n:
986        _init()
987        return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
988    if m == n + 1 or n == 0:
989        if m.mod(2).is_zero():
990            denom = (1 - x**2)**(m/2)
991        else:
992            denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)   
993        if m == n + 1:
994            return (-1)**m*(m-1).factorial()*2**n/denom
995        else:
996            return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
997    else:
998        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)
999
1000def hermite(n,x):
1001    """
1002    Returns the Hermite polynomial for integers `n > -1`.
1003   
1004    REFERENCE:
1005
1006    - AS 22.5.40 and 22.5.41, page 779.
1007   
1008    EXAMPLES::
1009   
1010        sage: x = PolynomialRing(QQ, 'x').gen()
1011        sage: hermite(2,x)
1012        4*x^2 - 2
1013        sage: hermite(3,x)
1014        8*x^3 - 12*x
1015        sage: hermite(3,2)
1016        40
1017        sage: S.<y> = PolynomialRing(RR)
1018        sage: hermite(3,y)
1019        8.00000000000000*y^3 - 12.0000000000000*y
1020        sage: R.<x,y> = QQ[]
1021        sage: hermite(3,y^2)
1022        8*y^6 - 12*y^2
1023        sage: w = var('w')
1024        sage: hermite(3,2*w)
1025        8*(8*w^2 - 3)*w
1026    """
1027    _init()
1028    return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
1029
1030def jacobi_P(n,a,b,x):
1031    r"""
1032    Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for
1033    integers `n > -1` and a and b symbolic or `a > -1`
1034    and `b > -1`. The Jacobi polynomials are actually defined
1035    for all a and b. However, the Jacobi polynomial weight
1036    `(1-x)^a(1+x)^b` isn't integrable for `a \leq -1`
1037    or `b \leq -1`.
1038   
1039    REFERENCE:
1040
1041    - table on page 789 in AS.
1042   
1043    EXAMPLES::
1044   
1045        sage: x = PolynomialRing(QQ, 'x').gen()
1046        sage: jacobi_P(2,0,0,x)
1047        3/2*x^2 - 1/2
1048        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
1049        5.009999999999998
1050    """
1051    _init()
1052    return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
1053
1054def laguerre(n,x):
1055    """
1056    Returns the Laguerre polynomial for integers `n > -1`.
1057   
1058    REFERENCE:
1059
1060    - AS 22.5.16, page 778 and AS page 789.
1061   
1062    EXAMPLES::
1063   
1064        sage: x = PolynomialRing(QQ, 'x').gen()
1065        sage: laguerre(2,x)
1066        1/2*x^2 - 2*x + 1
1067        sage: laguerre(3,x)
1068        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1069        sage: laguerre(2,2)
1070        -1
1071    """
1072    _init()
1073    return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
1074
1075
1076def legendre_Q(n,x):
1077    """
1078    Returns the Legendre function of the second kind for integers
1079    `n>-1`.
1080   
1081    Computed using Maxima.
1082   
1083    EXAMPLES::
1084   
1085        sage: P.<t> = QQ[]
1086        sage: legendre_Q(2, t)
1087        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
1088        sage: legendre_Q(3, 0.5)
1089        -0.198654771479482
1090        sage: legendre_Q(4, 2)
1091        443/16*I*pi + 443/16*log(3) - 365/12
1092        sage: legendre_Q(4, 2.0)
1093        0.00116107583162324 + 86.9828465962674*I
1094    """
1095    _init()
1096    return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
1097
1098def ultraspherical(n,a,x):
1099    """
1100    Returns the ultraspherical (or Gegenbauer) polynomial for integers
1101    `n > -1`.
1102   
1103    Computed using Maxima.
1104   
1105    REFERENCE:
1106
1107    - AS 22.5.27
1108   
1109    EXAMPLES::
1110   
1111        sage: x = PolynomialRing(QQ, 'x').gen()
1112        sage: ultraspherical(2,3/2,x)
1113        15/2*x^2 - 3/2
1114        sage: ultraspherical(2,1/2,x)
1115        3/2*x^2 - 1/2
1116        sage: ultraspherical(1,1,x)
1117        2*x     
1118        sage: t = PolynomialRing(RationalField(),"t").gen()
1119        sage: gegenbauer(3,2,t)
1120        32*t^3 - 12*t
1121    """
1122    _init()
1123    return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
1124
1125gegenbauer = ultraspherical