# Ticket #9706: orthogonal_polys.6.py

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

Version from 19. August 2010

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