Ticket #9706: orthogonal_polys.9.py

File orthogonal_polys.9.py, 71.8 KB (added by maldun, 11 years ago)

ortho polys with scipy support

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