Ticket #9706: orthogonal_polys.7.py

File orthogonal_polys.7.py, 79.3 KB (added by maldun, 11 years ago)

Added numpy support, eliminated some bugs (19.08.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#
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
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
379    maxima.eval('load("orthopoly");')
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
385#load /home/maldun/sage/sage-4.5.2/devel/sage-ortho/sage/functions/orthogonal_polys.py
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    #def _eval_numpy_(self, *args):
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        #Now this only evaluates the array pointwise, and only the first one...
570        #if isinstance(arg[0], numpy.ndarray)
571       
572       
573       
574       
575class Func_chebyshev_T(OrthogonalPolynomial): 
576
577    """
578    Class for the Chebyshev polynomial of the first kind.
579   
580    REFERENCE:
581
582    - AS 22.5.31 page 778 and AS 6.1.22 page 256.
583   
584    EXAMPLES::
585   
586       sage: chebyshev_T(3,x)
587       4*x^3 - 3*x
588       sage: var('k')
589       k
590       sage: test = chebyshev_T(k,x)
591       sage: test
592       chebyshev_T(k, x)
593    """
594   
595    def __init__(self):
596        OrthogonalPolynomial.__init__(self,"chebyshev_T",nargs = 2,
597conversions =dict(maxima='chebyshev_t',mathematica='ChebyshevT'))
598   
599    def _eval_special_values_(self,*args):
600        """
601        Values known for special values of x.
602        For details see A.S. 22.4 (p. 777)
603       
604        EXAMPLES:
605           
606            sage: var('n')
607            n
608            sage: chebyshev_T(n,1)
609            1
610            sage: chebyshev_T(n,-1)
611            (-1)^n
612        """
613        if args[-1] == 1:
614            return 1
615       
616        if args[-1] == -1:
617            return (-1)**args[0]
618
619        if (args[-1] == 0):
620            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
621
622        raise ValueError("Value not found!")
623   
624    def _evalf_(self, *args,**kwds):
625        """
626        Evals chebyshev_T
627        numerically with mpmath.
628        EXAMPLES::
629            sage: chebyshev_T(10,3).n(75)
630            2.261953700000000000000e7
631        """
632        try:
633            step_parent = kwds['parent']
634        except KeyError:
635            step_parent = parent(args[-1])
636
637        import sage.libs.mpmath.all as mpmath
638
639        try:
640            precision = step_parent.prec()
641        except AttributeError:
642            precision = mpmath.mp.prec
643
644        return mpmath.call(mpmath.chebyt,args[0],args[-1],prec = precision)
645
646    def _maxima_init_evaled_(self, *args):
647        n = args[0]
648        x = args[1]
649        return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
650
651    def _clenshaw_method_(self,*args):
652            """
653            Clenshaw method for chebyshev_T (means use recursions in this case)
654            This is much faster for numerical evaluation than maxima!
655            See A.S. 227 (p. 782) for details for the recurions
656            """
657
658            k = args[0]
659            x = args[1]
660
661            if k == 0:
662                return 1
663            elif k == 1:
664                return x
665            else:
666                #TODO: When evaluation of Symbolic Expressions works better
667                #use these explicit formulas instead!
668                #if -1 <= x <= 1:
669                #    return cos(k*acos(x))
670                #elif 1 < x:
671                #    return cosh(k*acosh(x))
672                #else: # x < -1
673                #    return (-1)**(k%2)*cosh(k*acosh(-x))
674               
675                help1 = 1
676                help2 = x
677                if is_Expression(x):
678                    #raise NotImplementedError
679                    help3 = 0
680                    for j in xrange(0,floor(k/2)+1):
681                        help3 = help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j-1)/factorial(j)/factorial(k-2*j)
682                    help3 = help3*k/2
683                else:
684                     for j in xrange(0,k-1):
685                        help3 = 2*x*help2 - help1
686                        help1 = help2
687                        help2 = help3
688               
689                return help3
690
691    def _eval_numpy_(self, *args):
692        #TODO: numpy_eval with help of the new scipy version!!!!
693        #Reason scipy 0.8 supports stable and fast numerical evaluation
694        #of ortho polys
695        #Now this only evaluates the array pointwise, and only the first one...
696        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
697        """
698        EXAMPLES::
699        sage: import numpy
700        sage: z = numpy.array([1,2])
701        sage: z2 = numpy.array([[1,2],[1,2]])
702        sage: z3 = numpy.array([1,2,3.])
703        sage: chebyshev_T(1,z)
704        array([1, 2])
705        sage: chebyshev_T(1,z2)
706        array([[1, 2],
707              [1, 2]])
708        sage: chebyshev_T(1,z3)
709        array([ 1.,  2.,  3.])
710
711        """
712        if isinstance(args[0],numpy.ndarray):
713            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
714
715        result = numpy.zeros(args[-1].shape).tolist()
716        if isinstance(args[-1][0],int):
717            for k in xrange(len(args[-1])):
718                result[k] = chebyshev_T(args[0],ZZ(args[-1][k]))
719
720        if isinstance(args[-1][0],float):
721            for k in xrange(len(args[-1])):
722                result[k] = chebyshev_T(args[0],RR(args[-1][k]))
723
724        if isinstance(args[-1][0],numpy.ndarray):
725            for k in xrange(len(args[-1])):
726                result[k] = chebyshev_T(args[0],args[-1][k])
727
728        return numpy.array(result)
729
730    def _derivative_(self, *args, **kwds):
731        """
732        Returns the derivative of chebyshev_T in form of the chebyshev Polynomial
733        of the second kind chebyshev_U
734        EXAMPLES::
735            sage: var('k')
736            k
737            sage: derivative(chebyshev_T(k,x),x)
738            k*chebyshev_U(k - 1, x)
739            sage: derivative(chebyshev_T(3,x),x)
740            12*x^2 - 3
741            sage: derivative(chebyshev_T(k,x),k)
742            Traceback (most recent call last):
743            ...
744            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
745 
746        """
747        diff_param = kwds['diff_param']
748        if diff_param == 0: 
749            raise NotImplementedError(
750"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
751        else:
752            return args[0]*chebyshev_U(args[0]-1,args[1])
753
754   
755chebyshev_T = Func_chebyshev_T()
756
757class Func_chebyshev_U(OrthogonalPolynomial): 
758   
759    """
760    Class for the Chebyshev polynomial of the second kind.
761   
762    REFERENCE:
763
764    - AS, 22.8.3 page 783 and AS 6.1.22 page 256.
765   
766    EXAMPLES::
767   
768        sage: chebyshev_U(3,x)
769        8*x^3 - 4*x
770    """
771
772    def __init__(self):
773        OrthogonalPolynomial.__init__(self,"chebyshev_U",
774nargs = 2,conversions =dict(maxima='chebyshev_u',mathematica='ChebyshevU'))
775       
776    def _clenshaw_method_(self,*args):
777        """
778        Clenshaw method for chebyshev_U (means use the recursion...)
779        This is much faster for numerical evaluation than maxima!
780        See A.S. 227 (p. 782) for details for the recurions
781        """
782        k = args[0]
783        x = args[1]
784           
785        if k == 0:
786            return 1
787        elif k == 1:
788            return 2*x
789        else:
790            help1 = 1
791            help2 = 2*x
792            if is_Expression(x):
793                #raise NotImplementedError("Maxima is faster here!")
794                help3 = 0
795                for j in xrange(0,floor(k/2)+1):
796                    help3 = help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j)/factorial(j)/factorial(k-2*j)
797               
798            else:
799                for j in xrange(0,k-1):
800                    help3 = 2*x*help2 - help1
801                    help1 = help2
802                    help2 = help3
803                   
804            return help3   
805   
806    def _maxima_init_evaled_(self, *args):
807        """
808        Uses 
809        EXAMPLES::
810        sage: x = PolynomialRing(QQ, 'x').gen()
811        sage: chebyshev_T(2,x)
812        2*x^2 - 1
813        """
814        n = args[0]
815        x = args[1]
816        return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
817
818       
819    def _evalf_(self, *args,**kwds):
820        """
821        Evals chebyshev_U
822        numerically with mpmath.
823        EXAMPLES::
824            sage: chebyshev_U(10,3).n(75)
825            4.661117900000000000000e7
826        """
827        try:
828            step_parent = kwds['parent']
829        except KeyError:
830            step_parent = parent(args[-1])
831
832        import sage.libs.mpmath.all as mpmath
833
834        try:
835            precision = step_parent.prec()
836        except AttributeError:
837            precision = mpmath.mp.prec
838
839        return mpmath.call(mpmath.chebyu,args[0],args[-1],prec = precision)
840
841    def _eval_special_values_(self,*args):
842        """
843        Special values known. A.S. 22.4 (p.777).
844        EXAMPLES:
845       
846            sage: var('n')
847            n
848            sage: chebyshev_U(n,1)
849            n + 1
850            sage: chebyshev_U(n,-1)
851            (n + 1)*(-1)^n
852        """
853        if args[-1] == 1:
854            return (args[0]+1)
855       
856        if args[-1] == -1:
857            return (-1)**args[0]*(args[0]+1)
858
859        if (args[-1] == 0):
860            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
861
862        raise ValueError("Value not found")
863
864    def _eval_numpy_(self, *args):
865        #TODO: numpy_eval with help of the new scipy version!!!!
866        #Reason scipy 0.8 supports stable and fast numerical evaluation
867        #of ortho polys
868        #Now this only evaluates the array pointwise, and only the first one...
869        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
870        """
871        EXAMPLES::
872        sage: import numpy
873        sage: z = numpy.array([1,2])
874        sage: z2 = numpy.array([[1,2],[1,2]])
875        sage: z3 = numpy.array([1,2,3.])
876        sage: chebyshev_U(1,z)
877        array([2, 4])
878        sage: chebyshev_U(1,z2)
879        array([[2, 4],
880              [2, 4]])
881        sage: chebyshev_U(1,z3)
882        array([ 2.,  4.,  6.])
883
884        """
885        if isinstance(args[0],numpy.ndarray):
886            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
887
888        result = numpy.zeros(args[-1].shape).tolist()
889        if isinstance(args[-1][0],int):
890            for k in xrange(len(args[-1])):
891                result[k] = chebyshev_U(args[0],ZZ(args[-1][k]))
892
893        if isinstance(args[-1][0],float):
894            for k in xrange(len(args[-1])):
895                result[k] = chebyshev_U(args[0],RR(args[-1][k]))
896
897        if isinstance(args[-1][0],numpy.ndarray):
898            for k in xrange(len(args[-1])):
899                result[k] = chebyshev_U(args[0],args[-1][k])
900
901        return numpy.array(result)
902
903
904    def _derivative_(self, *args, **kwds):
905        """
906        Returns the derivative of chebyshev_U in form of the chebyshev
907        Polynomials of the first and second kind
908       
909        EXAMPLES::
910            sage: var('k')
911            k
912            sage: derivative(chebyshev_U(k,x),x)
913            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
914            sage: derivative(chebyshev_U(3,x),x)
915            24*x^2 - 4
916            sage: derivative(chebyshev_U(k,x),k)
917            Traceback (most recent call last):
918            ...
919            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
920 
921        """
922        diff_param = kwds['diff_param']
923        if diff_param == 0: 
924            raise NotImplementedError(
925"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
926        else:
927            return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
928                    chebyshev_U(args[0],args[1]))/(args[1]**2-1)
929       
930chebyshev_U = Func_chebyshev_U()
931
932class Func_gegenbauer(OrthogonalPolynomial):
933    """
934    Returns the ultraspherical (or Gegenbauer) polynomial.
935   
936    REFERENCE:
937
938    - AS 22.5.27
939   
940    EXAMPLES::
941   
942        sage: x = PolynomialRing(QQ, 'x').gen()
943        sage: ultraspherical(2,3/2,x)
944        15/2*x^2 - 3/2
945        sage: ultraspherical(2,1/2,x)
946        3/2*x^2 - 1/2
947        sage: ultraspherical(1,1,x)
948        2*x     
949        sage: t = PolynomialRing(RationalField(),"t").gen()
950        sage: gegenbauer(3,2,t)
951        32*t^3 - 12*t
952    """   
953    def __init__(self):
954        OrthogonalPolynomial.__init__(self,"gegenbauer",nargs = 3,conversions =dict(maxima='ultraspherical',mathematica='GegenbauerC'))
955
956    def _clenshaw_method_(self,*args):
957        """
958        Clenshaw method for gegenbauer poly (means use the recursion...)
959        This is much faster for numerical evaluation than maxima!
960        See A.S. 227 (p. 782) for details for the recurions
961        """
962        k = args[0]
963        x = args[-1]
964        alpha = args[1]
965
966        if is_Expression(alpha) or abs(alpha) > numpy.finfo(float).eps:
967            alpha_zero = False
968        else:
969            alpha_zero = True
970
971        if k == 0:
972            return 1
973        elif k == 1:
974            if not alpha_zero:
975                return 2*x*alpha
976            else:
977                return 2*x # It seems that maxima evals this the wrong way! (see A.S. 22.4 (p.777))
978        else:
979            help1 = 1
980            if alpha_zero:
981                help2 = 2*x
982                gamma_alpha = 1
983            else:
984                help2 = 2*x*alpha
985                gamma_alpha = gamma(alpha)
986
987            help3 = 0
988            if is_Expression(x):
989                for j in xrange(0,floor(k/2)+1):
990                    help3 = help3 + (-1)**j*gamma(alpha+k-j)/factorial(j)/factorial(k-2*j)*(2*x)**(k-2*j)/gamma_alpha
991               
992            else:
993                for j in xrange(1,k):
994                    help3 = 2*(j+alpha)*x*help2 - (j+2*alpha-1)*help1
995                    help3 = help3/(j+1)
996                    help1 = help2
997                    help2 = help3
998                   
999            return help3   
1000
1001    def _maxima_init_evaled_(self, *args):
1002        """
1003        Uses 
1004        EXAMPLES::
1005        sage: x = PolynomialRing(QQ, 'x').gen()
1006        sage: ultraspherical(2,1/2,x)
1007        3/2*x^2 - 1/2
1008        """
1009        n = args[0]
1010        a = args[1]
1011        x = args[2]
1012        return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
1013
1014
1015    def _evalf_(self, *args,**kwds):
1016        """
1017        Evals gegenbauer
1018        numerically with mpmath.
1019        EXAMPLES::
1020            sage: gegenbauer(10,2,3.).n(54)
1021            5.25360702000000e8
1022        """
1023        if args[1] == 0:
1024            raise ValueError("mpmath don't handle alpha = 0 correctly!")
1025           
1026        try:
1027            step_parent = kwds['parent']
1028        except KeyError:
1029            step_parent = parent(args[-1])
1030
1031        import sage.libs.mpmath.all as mpmath
1032
1033        try:
1034            precision = step_parent.prec()
1035        except AttributeError:
1036            precision = mpmath.mp.prec
1037
1038        return mpmath.call(mpmath.gegenbauer,args[0],args[1],args[-1],prec = precision)
1039
1040    def _eval_special_values_(self,*args):
1041        """
1042        Special values known. A.S. 22.4 (p.777)
1043        EXAMPLES:
1044       
1045            sage: var('n a')
1046            (n, a)
1047            sage: gegenbauer(n,1/2,x)
1048            legendre_P(n, x)
1049            sage: gegenbauer(n,0,x)
1050            1/2*n*chebyshev_T(n, x)
1051            sage: gegenbauer(n,1,x)
1052            chebyshev_U(n, x)
1053            sage: gegenbauer(n,a,1)
1054            binomial(2*a + n - 1, n)
1055            sage: gegenbauer(n,a,-1)
1056            (-1)^n*binomial(2*a + n - 1, n)
1057            sage: gegenbauer(n,a,0)
1058            1/2*((-1)^n + 1)*(-1)^(1/2*n)*gamma(a + 1/2*n)/(gamma(a)*gamma(1/2*n))
1059        """
1060        if args[1] == 0 and args[0] != 0:
1061            return args[0]*chebyshev_T(args[0],args[-1])/2
1062
1063        if args[1] == 0.5:
1064            return legendre_P(args[0],args[-1])
1065
1066        if args[1] == 1:
1067            return chebyshev_U(args[0],args[-1])
1068
1069        if args[-1] == 1:
1070            return binomial(args[0] + 2*args[1] - 1,args[0])
1071       
1072        if args[-1] == -1:
1073            return (-1)**args[0]*binomial(args[0] + 2*args[1] - 1,args[0])
1074
1075        if (args[-1] == 0):
1076            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2*gamma(args[1]+args[0]/2)/gamma(args[1])/gamma(args[0]/2)
1077
1078        raise ValueError("Value not found")
1079
1080    def _eval_numpy_(self, *args):
1081        #TODO: numpy_eval with help of the new scipy version!!!!
1082        #Reason scipy 0.8 supports stable and fast numerical evaluation
1083        #of ortho polys
1084        #Now this only evaluates the array pointwise, and only the first one...
1085        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
1086        """
1087        EXAMPLES::
1088        sage: import numpy
1089        sage: z = numpy.array([1,2])
1090        sage: z2 = numpy.array([[1,2],[1,2]])
1091        sage: z3 = numpy.array([1,2,3.])
1092        sage: gegenbauer(2,0,z)
1093        array([1, 7])
1094        sage: gegenbauer(2,1,z)
1095        array([ 3, 15])
1096        sage: gegenbauer(2,1,z2)
1097        array([[ 3, 15],
1098              [ 3, 15]])
1099        sage: gegenbauer(2,1,z3)
1100        array([  3.,  15.,  35.])
1101        """
1102
1103        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
1104            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
1105       
1106        result = numpy.zeros(args[-1].shape).tolist()
1107        if isinstance(args[-1][0],int):
1108            for k in xrange(len(args[-1])):
1109                result[k] = gegenbauer(args[0],args[1],ZZ(args[-1][k]))
1110
1111        if isinstance(args[-1][0],float):
1112            for k in xrange(len(args[-1])):
1113                result[k] = gegenbauer(args[0],args[1],RR(args[-1][k]))
1114
1115        if isinstance(args[-1][0],numpy.ndarray):
1116            for k in xrange(len(args[-1])):
1117                result[k] = gegenbauer(args[0],args[1],args[-1][k])
1118
1119        return numpy.array(result)
1120
1121    def _derivative_(self, *args, **kwds):
1122        """
1123        Returns the derivative of chebyshev_U in form of the chebyshev
1124        Polynomials of the first and second kind
1125       
1126        EXAMPLES::
1127            sage: var('k a')
1128            (k, a)
1129            sage: derivative(gegenbauer(k,a,x),x)
1130            (k*x*gegenbauer(k, a, x) - (2*a + k - 1)*gegenbauer(k - 1, a, x))/(x^2 - 1)
1131            sage: derivative(gegenbauer(4,3,x),x)
1132            960*x^3 - 240*x
1133            sage: derivative(gegenbauer(k,a,x),a)
1134            Traceback (most recent call last):
1135            ...
1136            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
1137 
1138        """
1139        diff_param = kwds['diff_param']
1140        if diff_param in [0,1]: 
1141            raise NotImplementedError(
1142"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1143        else:
1144            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)
1145
1146gegenbauer = Func_gegenbauer()
1147ultraspherical = Func_gegenbauer()
1148
1149class Func_gen_laguerre(OrthogonalPolynomial):
1150    """
1151    Returns the generalized Laguerre polynomial for integers `n > -1`.
1152    Typically, a = 1/2 or a = -1/2.
1153   
1154    REFERENCE:
1155
1156    - table on page 789 in AS.
1157   
1158    EXAMPLES::
1159   
1160        sage: x = PolynomialRing(QQ, 'x').gen()
1161        sage: gen_laguerre(2,1,x)
1162        1/2*x^2 - 3*x + 3
1163        sage: gen_laguerre(2,1/2,x)
1164        1/2*x^2 - 5/2*x + 15/8
1165        sage: gen_laguerre(2,-1/2,x)
1166        1/2*x^2 - 3/2*x + 3/8
1167        sage: gen_laguerre(2,0,x)
1168        1/2*x^2 - 2*x + 1
1169        sage: gen_laguerre(3,0,x)
1170        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1171    """
1172    def __init__(self):
1173        OrthogonalPolynomial.__init__(self,"gen_laguerre",nargs = 3,
1174conversions =dict(maxima='gen_laguerre',mathematica='LaguerreL'))
1175       
1176    def _clenshaw_method_(self,*args):
1177        """
1178        Clenshaw method for gen_laguerre polynomial (means use the recursion...)
1179        This is much faster for numerical evaluation than maxima!
1180        See A.S. 227 (p. 782) for details for the recurions.
1181        Warning: The clanshaw method for the laguerre Polynomials
1182        should only used for exact data types, when high orders are
1183        used, due to weak instabilities of the recursion!
1184        """
1185        k = args[0]
1186        a = args[1]
1187        x = args[-1]
1188       
1189        if k == 0:
1190            return 1
1191        elif k == 1:
1192            return 1-x + a
1193        else:
1194            help1 = 1
1195            help2 = 1-x + a
1196            if is_Expression(x): 
1197                help3 = 0
1198                for j in xrange(0,k+1):
1199                    help3 = help3 + (-1)**j*binomial(k+a,k-j)/factorial(j)*x**j
1200            else:
1201                for j in xrange(1,k):
1202                    help3 = -x*help2 + (2*j+1+a)*help2 - (j+a)*help1
1203                    help3 = help3/(j+1)
1204                    help1 = help2
1205                    help2 = help3
1206               
1207            return help3   
1208
1209    def _maxima_init_evaled_(self, *args):
1210        n = args[0]
1211        a = args[1]
1212        x = args[-1]
1213       
1214        from sage.functions.all import sqrt
1215        _init()
1216        return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
1217
1218    def _evalf_(self, *args,**kwds):
1219        """
1220        Evals laguerre polynomial
1221        numerically with mpmath.
1222        EXAMPLES::
1223            sage: laguerre(3,5.).n(53)
1224            2.66666666666667
1225        """
1226        try:
1227            step_parent = kwds['parent']
1228        except KeyError:
1229            step_parent = parent(args[-1])
1230
1231        import sage.libs.mpmath.all as mpmath
1232
1233        try:
1234            precision = step_parent.prec()
1235        except AttributeError:
1236            precision = mpmath.mp.prec
1237
1238        return mpmath.call(mpmath.laguerre,args[0],args[1],args[-1],prec = precision)
1239   
1240    def _eval_special_values_(self,*args):
1241        """
1242        Special values known.
1243        EXAMPLES:
1244           
1245            sage: var('n')
1246            n
1247            sage: gen_laguerre(n,1,0)
1248            binomial(n + 1, n)
1249        """
1250       
1251        if (args[-1] == 0):
1252            try:
1253                return binomial(args[0]+args[1],args[0])
1254            except TypeError:
1255                pass
1256
1257        raise ValueError("Value not found")
1258
1259    def _eval_numpy_(self, *args):
1260        #TODO: numpy_eval with help of the new scipy version!!!!
1261        #Reason scipy 0.8 supports stable and fast numerical evaluation
1262        #of ortho polys
1263        #Now this only evaluates the array pointwise, and only the first one...
1264        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
1265        """
1266        EXAMPLES::
1267        sage: import numpy
1268        sage: z = numpy.array([1,2])
1269        sage: z2 = numpy.array([[1,2],[1,2]])
1270        sage: z3 = numpy.array([1,2,3.])
1271        sage: gen_laguerre(1,1,z)
1272        array([1, 0])
1273        sage: gen_laguerre(1,1,z2)
1274        array([[1, 0],
1275              [1, 0]])
1276        sage: gen_laguerre(1,1,z3)
1277        array([ 1.,  0., -1.])
1278        """
1279
1280        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
1281            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
1282       
1283        result = numpy.zeros(args[-1].shape).tolist()
1284        if isinstance(args[-1][0],int):
1285            for k in xrange(len(args[-1])):
1286                result[k] = gen_laguerre(args[0],args[1],ZZ(args[-1][k]))
1287
1288        if isinstance(args[-1][0],float):
1289            for k in xrange(len(args[-1])):
1290                result[k] = gen_laguerre(args[0],args[1],RR(args[-1][k]))
1291
1292        if isinstance(args[-1][0],numpy.ndarray):
1293            for k in xrange(len(args[-1])):
1294                result[k] = gen_laguerre(args[0],args[1],args[-1][k])
1295
1296        return numpy.array(result)
1297
1298   
1299
1300    def _derivative_(self,*args,**kwds):
1301        """return the derivative of gen_laguerre in
1302        form of the gen. Laguerre Polynomial.
1303        EXAMPLES::
1304        sage: n = var('n')
1305        sage: derivative(gen_laguerre(3,1,x),x)
1306        -1/2*x^2 + 4*x - 6
1307        sage: derivative(gen_laguerre(n,1,x),x)
1308        -((n + 1)*gen_laguerre(n - 1, 1, x) - n*gen_laguerre(n, 1, x))/x
1309        """
1310        diff_param = kwds['diff_param']
1311        if diff_param == 0: 
1312            raise NotImplementedError(
1313"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1314        else:
1315            return (args[0]*gen_laguerre(args[0],args[1],args[-1])-(args[0]+args[1])*gen_laguerre(args[0]-1,args[1],args[-1]))/args[-1]
1316
1317gen_laguerre = Func_gen_laguerre()
1318
1319class Func_hermite(OrthogonalPolynomial):
1320    """
1321    Returns the Hermite polynomial for integers `n > -1`.
1322   
1323    REFERENCE:
1324
1325    - AS 22.5.40 and 22.5.41, page 779.
1326   
1327    #sage: S.<y> = PolynomialRing(RR)<- Here a bug of the base class BuiltinFunction occours!
1328
1329    EXAMPLES::
1330   
1331        sage: x = PolynomialRing(QQ, 'x').gen()
1332        sage: hermite(2,x)
1333        4*x^2 - 2
1334        sage: hermite(3,x)
1335        8*x^3 - 12*x
1336        sage: hermite(3,2)
1337        40
1338        sage: y = var('y')
1339        sage: hermite(3,y).polynomial(PolynomialRing(RR, 'y'))
1340        8.00000000000000*y^3 - 12.0000000000000*y
1341        sage: R.<x,y> = QQ[]
1342        sage: hermite(3,y^2)
1343        8*y^6 - 12*y^2
1344        sage: w = var('w')
1345        sage: hermite(3,2*w)
1346        64*w^3 - 24*w
1347    """
1348   
1349    def __init__(self):
1350        OrthogonalPolynomial.__init__(self,"hermite",nargs = 2,
1351conversions =dict(maxima='hermite',mathematica='HermiteH'))
1352
1353    def _clenshaw_method_(self,*args):
1354        """
1355        Clenshaw method for hermite polynomial (means use the recursion...)
1356        See A.S. 227 (p. 782) for details for the recurions
1357        For the symbolic evaluation, maxima seems to be quite fast.
1358        The break even point between the recursion and Maxima is about
1359        n = 25
1360        """
1361        k = args[0]
1362        x = args[1]
1363           
1364        if k == 0:
1365            return 1
1366        elif k == 1:
1367            return 2*x
1368        else:
1369            help1 = 1
1370            help2 = 2*x
1371            if is_Expression(x):
1372                help3 = 0
1373                for j in xrange(0,floor(k/2)+1):
1374                    help3 = help3 + (-1)**j*(2*x)**(k-2*j)/factorial(j)/factorial(k-2*j)
1375                help3 = help3*factorial(k)
1376            else:
1377                for j in xrange(1,k):
1378                    help3 = 2*x*help2 - 2*j*help1
1379                    help1 = help2
1380                    help2 = help3
1381                   
1382            return help3   
1383   
1384           
1385    def _evalf_(self, *args,**kwds):
1386        """
1387        Evals hermite
1388        numerically with mpmath.
1389        EXAMPLES::
1390            sage: hermite(3,2.).n(74)
1391            40.0000000000000000000
1392        """
1393        try:
1394            step_parent = kwds['parent']
1395        except KeyError:
1396            step_parent = parent(args[-1])
1397
1398        import sage.libs.mpmath.all as mpmath
1399
1400        try:
1401            precision = step_parent.prec()
1402        except AttributeError:
1403            precision = mpmath.mp.prec
1404
1405        return mpmath.call(mpmath.hermite,args[0],args[-1],prec = precision)
1406   
1407    def _eval_special_values_(self,*args):
1408        """
1409        Special values known. A.S. 22.4 (p.777)
1410        EXAMPLES:
1411       
1412            sage: var('n')
1413            n
1414            sage: hermite(n,0)
1415            ((-1)^n + 1)*(-1)^(1/2*n)*factorial(n)/gamma(1/2*n + 1)
1416        """
1417
1418        if (args[-1] == 0):
1419            return (1+(-1)**args[0])*(-1)**(args[0]/2)*factorial(args[0])/gamma(args[0]/2+1)
1420
1421        raise ValueError("Value not found")
1422
1423    def _maxima_init_evaled_(self, *args):
1424        """
1425        Old maxima method.
1426        """
1427        n = args[0]
1428        x = args[1]
1429        return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
1430
1431    def _eval_numpy_(self, *args):
1432        #TODO: numpy_eval with help of the new scipy version!!!!
1433        #Reason scipy 0.8 supports stable and fast numerical evaluation
1434        #of ortho polys
1435        #Now this only evaluates the array pointwise, and only the first one...
1436        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
1437        """
1438        EXAMPLES::
1439        sage: import numpy
1440        sage: z = numpy.array([1,2])
1441        sage: z2 = numpy.array([[1,2],[1,2]])
1442        sage: z3 = numpy.array([1,2,3.])
1443        sage: hermite(1,z)
1444        array([2, 4])
1445        sage: hermite(1,z2)
1446        array([[2, 4],
1447              [2, 4]])
1448        sage: hermite(1,z3)
1449        array([ 2.,  4.,  6.])
1450
1451        """
1452        if isinstance(args[0],numpy.ndarray):
1453            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
1454
1455        result = numpy.zeros(args[-1].shape).tolist()
1456        if isinstance(args[-1][0],int):
1457            for k in xrange(len(args[-1])):
1458                result[k] = hermite(args[0],ZZ(args[-1][k]))
1459
1460        if isinstance(args[-1][0],float):
1461            for k in xrange(len(args[-1])):
1462                result[k] = hermite(args[0],RR(args[-1][k]))
1463
1464        if isinstance(args[-1][0],numpy.ndarray):
1465            for k in xrange(len(args[-1])):
1466                result[k] = hermite(args[0],args[-1][k])
1467
1468        return numpy.array(result)
1469
1470
1471    def _derivative_(self, *args, **kwds):
1472        """
1473        Returns the derivative of the hermite polynomial in form of the chebyshev
1474        Polynomials of the first and second kind
1475       
1476        EXAMPLES::
1477            sage: var('k')
1478            k
1479            sage: derivative(hermite(k,x),x)
1480            2*k*hermite(k - 1, x)
1481             
1482        """
1483        diff_param = kwds['diff_param']
1484        if diff_param == 0: 
1485            raise NotImplementedError(
1486"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1487        else:
1488            return 2*args[0]*hermite(args[0]-1,args[1])
1489 
1490hermite = Func_hermite()       
1491
1492
1493class Func_jacobi_P(OrthogonalPolynomial):
1494    r"""
1495    Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for
1496    integers `n > -1` and a and b symbolic or `a > -1`
1497    and `b > -1`. The Jacobi polynomials are actually defined
1498    for all a and b. However, the Jacobi polynomial weight
1499    `(1-x)^a(1+x)^b` isn't integrable for `a \leq -1`
1500    or `b \leq -1`.
1501   
1502    REFERENCE:
1503
1504    - table on page 789 in AS.
1505   
1506    EXAMPLES::
1507   
1508        sage: x = PolynomialRing(QQ, 'x').gen()
1509        sage: jacobi_P(2,0,0,x)
1510        3/2*x^2 - 1/2
1511        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
1512        5.009999999999998
1513    """
1514
1515    def __init__(self):
1516        OrthogonalPolynomial.__init__(self,"jacobi_P",nargs = 4,
1517conversions =dict(maxima='jacobi_p',mathematica='JacobiP'))
1518   
1519    def _clenshaw_method_(self,*args):
1520        """
1521        Clenshaw method for jacobi_P (means use the recursion,
1522        or sum formula)
1523        This is much faster for numerical evaluation than maxima!
1524        See A.S. 227 (p. 782) for details for the recurions.
1525        Warning: The clanshaw method for the Jacobi Polynomials
1526        should only used for exact data types, when high orders are
1527        used, due to weak instabilities of the recursion!
1528        """ 
1529        k = args[0]
1530        x = args[-1]
1531        alpha = args[1]
1532        beta = args[2]
1533
1534        if k == 0:
1535            return 1
1536        elif k == 1:
1537            return (alpha-beta + (alpha+beta+2)*x)/2
1538        else:
1539           
1540            if is_Expression(x) or is_Expression(alpha) or is_Expression(beta):
1541                #Here we use the sum formula of jacobi_P it seems this is rather
1542                #optimal for use.
1543                help1 = gamma(alpha+k+1)/factorial(k)/gamma(alpha+beta+k+1)
1544                help2 = 0
1545                for j in xrange(0,k+1):
1546                    help2 = help2 + binomial(k,j)*gamma(alpha+beta+k+j+1)/gamma(alpha+j+1)*((x-1)/2)**j
1547                return help1*help2
1548            else:
1549                help1 = 1
1550                help2 = (alpha-beta + (alpha+beta+2)*x)/2
1551
1552                for j in xrange(1,k):
1553                    a1 = 2*(j+1)*(j+alpha+beta+1)*(2*j+alpha+beta)
1554                    a2 = (2*j+alpha+beta+1)*(alpha**2-beta**2)
1555                    a3 = (2*j+alpha+beta)*(2*j+alpha+beta+1)*(2*j+alpha+beta+2)
1556                    a4 = 2*(j+alpha)*(j+beta)*(2*j+alpha+beta+2)
1557                    help3 = (a2+a3*x)*help2 - a4*help1
1558                    help3 = help3/a1
1559                    help1 = help2
1560                    help2 = help3
1561               
1562            return help3   
1563
1564    def _maxima_init_evaled_(self, *args):
1565        """
1566        Old maxima method.
1567        """
1568        n = args[0]
1569        a = args[1]
1570        b = args[2]
1571        x = args[3]
1572        return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
1573
1574    def _evalf_(self, *args,**kwds):
1575        """
1576        Evals jacobi_P
1577        numerically with mpmath.
1578        EXAMPLES::
1579            sage: jacobi_P(10,2,3,3).n(75)
1580            1.322776620000000000000e8
1581        """
1582        try:
1583            step_parent = kwds['parent']
1584        except KeyError:
1585            step_parent = parent(args[-1])
1586
1587        import sage.libs.mpmath.all as mpmath
1588
1589        try:
1590            precision = step_parent.prec()
1591        except AttributeError:
1592            precision = mpmath.mp.prec
1593
1594        return mpmath.call(mpmath.jacobi,args[0],args[1],args[2],args[-1],prec = precision)
1595
1596    def _eval_special_values_(self,*args):
1597        """
1598        Special values known. A.S. 22.4 (p.777)
1599        EXAMPLES:
1600           
1601            sage: var('n')
1602            n
1603            sage: legendre_P(n,1)
1604            1
1605            sage: legendre_P(n,-1)
1606            (-1)^n
1607        """
1608        if args[-1] == 1:
1609            return binomial(args[0]+args[1],args[0])
1610       
1611        if args[-1] == -1:
1612            return (-1)**args[0]*binomial(args[0]+args[1],args[0])
1613       
1614        if args[1] == 0 and args[2] == 0:
1615            return legendre_P(args[0],args[-1])
1616
1617        if args[1] == -0.5 and args[2] == -0.5:
1618            try:
1619                return binomial(2*args[0],args[0])*chebyshev_T(args[0],args[-1])/4**args[0]
1620            except TypeError:
1621                pass
1622
1623        raise ValueError("Value not found")
1624
1625    def _eval_numpy_(self, *args):
1626        #TODO: numpy_eval with help of the new scipy version!!!!
1627        #Reason scipy 0.8 supports stable and fast numerical evaluation
1628        #of ortho polys
1629        #Now this only evaluates the array pointwise, and only the first one...
1630        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
1631        """
1632        EXAMPLES::
1633        sage: import numpy
1634        sage: z = numpy.array([1,2])
1635        sage: z2 = numpy.array([[1,2],[1,2]])
1636        sage: z3 = numpy.array([1,2,3.])
1637        sage: gen_laguerre(1,1,z)
1638        array([1, 0])
1639        sage: gen_laguerre(1,1,z2)
1640        array([[1, 0],
1641              [1, 0]])
1642        sage: gen_laguerre(1,1,z3)
1643        array([ 1.,  0., -1.])
1644        """
1645
1646        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray) or isinstance(args[2],numpy.ndarray):
1647            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
1648       
1649        result = numpy.zeros(args[-1].shape).tolist()
1650        if isinstance(args[-1][0],int):
1651            for k in xrange(len(args[-1])):
1652                result[k] = jacobi_P(args[0],args[1],args[2],ZZ(args[-1][k]))
1653
1654        if isinstance(args[-1][0],float):
1655            for k in xrange(len(args[-1])):
1656                result[k] = jacobi_P(args[0],args[1],args[2],RR(args[-1][k]))
1657
1658        if isinstance(args[-1][0],numpy.ndarray):
1659            for k in xrange(len(args[-1])):
1660                result[k] = jacobi_P(args[0],args[1],args[2],args[-1][k])
1661
1662        return numpy.array(result)
1663
1664    def _derivative_(self, *args, **kwds):
1665        """
1666        Returns the derivative of jacobi_P in form of jacobi_polynomials
1667       
1668        EXAMPLES::
1669            sage: var('k a b')
1670            (k, a, b)
1671            sage: derivative(jacobi_P(k,a,b,x),x)
1672            -(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))
1673            sage: derivative(jacobi_P(2,1,3,x),x)
1674            14*x - 7/2
1675            sage: derivative(jacobi_P(k,a,b,x),a)
1676            Traceback (most recent call last):
1677            ...
1678            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
1679 
1680        """
1681        diff_param = kwds['diff_param']
1682        if diff_param in [0,1,2]: 
1683            raise NotImplementedError(
1684"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1685        else:
1686            return (args[0]*(args[1]-args[2]-(2*args[0]+args[1]+args[2])*args[-1])*
1687                    jacobi_P(args[0],args[1],args[2],args[-1])+ 2*(args[0]+args[1])*
1688                    (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)
1689       
1690       
1691jacobi_P = Func_jacobi_P()       
1692
1693class Func_laguerre(OrthogonalPolynomial):
1694    """
1695    Returns the Laguerre polynomial for integers `n > -1`.
1696   
1697    REFERENCE:
1698
1699    - AS 22.5.16, page 778 and AS page 789.
1700   
1701    EXAMPLES::
1702   
1703        sage: x = PolynomialRing(QQ, 'x').gen()
1704        sage: laguerre(2,x)
1705        1/2*x^2 - 2*x + 1
1706        sage: laguerre(3,x)
1707        -1/6*x^3 + 3/2*x^2 - 3*x + 1
1708        sage: laguerre(2,2)
1709        -1
1710    """
1711    def __init__(self):
1712        OrthogonalPolynomial.__init__(self,"laguerre",nargs = 2,
1713conversions =dict(maxima='laguerre',mathematica='LaguerreL'))
1714       
1715    def _clenshaw_method_(self,*args):
1716        """
1717        Clenshaw method for laguerre polynomial (means use the recursion...)
1718        This is much faster for numerical evaluation than maxima!
1719        See A.S. 227 (p. 782) for details for the recurions.
1720        Warning: The clanshaw method for the laguerre Polynomials
1721        should only used for exact data types, when high orders are
1722        used, due to weak instabilities of the recursion!
1723        """
1724        k = args[0]
1725        x = args[-1]
1726       
1727        if k == 0:
1728            return 1
1729        elif k == 1:
1730            return 1-x
1731        else:
1732            help1 = 1
1733            help2 = 1-x
1734            if is_Expression(x):
1735                help3 = 0
1736                for j in xrange(0,k+1):
1737                    help3 = help3 + (-1)**j*binomial(k,k-j)/factorial(j)*x**j
1738            else:
1739                for j in xrange(1,k):
1740                    help3 = -x*help2 + (2*j+1)*help2 - j*help1
1741                    help3 = help3/(j+1)
1742                    help1 = help2
1743                    help2 = help3
1744               
1745            return help3   
1746
1747    def _maxima_init_evaled_(self, *args):
1748        n = args[0]
1749        x = args[1]
1750        return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
1751
1752    def _evalf_(self, *args,**kwds):
1753        """
1754        Evals laguerre polynomial
1755        numerically with mpmath.
1756        EXAMPLES::
1757            sage: laguerre(3,5.).n(53)
1758            2.66666666666667
1759        """
1760        try:
1761            step_parent = kwds['parent']
1762        except KeyError:
1763            step_parent = parent(args[-1])
1764
1765        import sage.libs.mpmath.all as mpmath
1766
1767        try:
1768            precision = step_parent.prec()
1769        except AttributeError:
1770            precision = mpmath.mp.prec
1771
1772        return mpmath.call(mpmath.laguerre,args[0],0,args[-1],prec = precision)
1773   
1774    def _eval_special_values_(self,*args):
1775        """
1776        Special values known.
1777        EXAMPLES:
1778           
1779            sage: var('n')
1780            n
1781            sage: laguerre(n,0)
1782            1
1783        """
1784       
1785        if (args[-1] == 0):
1786            try:
1787                return 1
1788            except TypeError:
1789                pass
1790
1791        raise ValueError("Value not found")
1792 
1793
1794    def _eval_numpy_(self, *args):
1795        #TODO: numpy_eval with help of the new scipy version!!!!
1796        #Reason scipy 0.8 supports stable and fast numerical evaluation
1797        #of ortho polys
1798        #Now this only evaluates the array pointwise, and only the first one...
1799        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
1800        """
1801        EXAMPLES::
1802        sage: import numpy
1803        sage: z = numpy.array([1,2])
1804        sage: z2 = numpy.array([[1,2],[1,2]])
1805        sage: z3 = numpy.array([1,2,3.])
1806        sage: laguerre(1,z)
1807        array([ 0, -1])
1808        sage: laguerre(1,z2)
1809        array([[ 0, -1],
1810               [ 0, -1]])
1811        sage: laguerre(1,z3)
1812        array([ 0., -1., -2.])
1813
1814        """
1815        if isinstance(args[0],numpy.ndarray):
1816            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
1817
1818        result = numpy.zeros(args[-1].shape).tolist()
1819        if isinstance(args[-1][0],int):
1820            for k in xrange(len(args[-1])):
1821                result[k] = laguerre(args[0],ZZ(args[-1][k]))
1822
1823        if isinstance(args[-1][0],float):
1824            for k in xrange(len(args[-1])):
1825                result[k] = laguerre(args[0],RR(args[-1][k]))
1826
1827        if isinstance(args[-1][0],numpy.ndarray):
1828            for k in xrange(len(args[-1])):
1829                result[k] = laguerre(args[0],args[-1][k])
1830
1831        return numpy.array(result)
1832
1833    def _derivative_(self,*args,**kwds):
1834        """return the derivative of laguerre in
1835        form of the Laguerre Polynomial.
1836        EXAMPLES::
1837        sage: n = var('n')
1838        sage: derivative(laguerre(3,x),x)
1839        -1/2*x^2 + 3*x - 3
1840        sage: derivative(laguerre(n,x),x)
1841        -(n*laguerre(n - 1, x) - n*laguerre(n, x))/x
1842        """
1843        diff_param = kwds['diff_param']
1844        if diff_param == 0: 
1845            raise NotImplementedError(
1846"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
1847        else:
1848            return (args[0]*laguerre(args[0],args[-1])-args[0]*laguerre(args[0]-1,args[-1]))/args[1]
1849
1850laguerre = Func_laguerre()
1851
1852class Func_legendre_P(OrthogonalPolynomial):
1853    """
1854    Returns the Legendre polynomial of the first kind..
1855   
1856    REFERENCE:
1857
1858    - AS 22.5.35 page 779.
1859   
1860    EXAMPLES::
1861   
1862        sage: P.<t> = QQ[]
1863        sage: legendre_P(2,t)
1864        3/2*t^2 - 1/2
1865        sage: legendre_P(3, 1.1)
1866        1.67750000000000
1867        sage: legendre_P(3, 1 + I)
1868        7/2*I - 13/2
1869        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
1870        Traceback (most recent call last):
1871        ...
1872        TypeError: cannot coerce arguments: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Integer Ring to Symbolic Ring
1873        sage: legendre_P(3, GF(11)(5))
1874        8
1875    """
1876    def __init__(self):
1877        OrthogonalPolynomial.__init__(self,"legendre_P",nargs = 2,
1878conversions =dict(maxima='legendre_p',mathematica='LegendreP'))
1879
1880    def _clenshaw_method_(self,*args):
1881        """
1882        Clenshaw method for legendre_P (means use the recursion...)
1883        This is much faster for numerical evaluation than maxima!
1884        See A.S. 227 (p. 782) for details for the recurions.
1885        Warning: The clanshaw method for the Legendre Polynomials
1886        should only used for exact data types, when high orders are
1887        used, due to weak instabilities of the recursion!
1888        """
1889        k = args[0]
1890        x = args[-1]
1891       
1892        if k == 0:
1893            return 1
1894        elif k == 1:
1895            return x
1896        else:
1897            help1 = 1
1898            help2 = x
1899            if is_Expression(x):
1900                #raise NotImplementedError("Maxima is faster here...")
1901                help1 = ZZ(2**k) #Workarround because of segmentation fault...
1902                help3 = 0
1903                for j in xrange(0,floor(k/2)+1):
1904                    help3 = help3 + (-1)**j*x**(k-2*j)*binomial(k,j)*binomial(2*(k-j),k)
1905               
1906                help3 = help3/help1
1907            else:
1908                for j in xrange(1,k):
1909                    help3 = (2*j+1)*x*help2 - j*help1
1910                    help3 = help3/(j+1)
1911                    help1 = help2
1912                    help2 = help3
1913               
1914            return help3   
1915
1916    def _maxima_init_evaled_(self, *args):
1917        n = args[0]
1918        x = args[1]
1919        return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)), locals={'x':x})
1920
1921    def _evalf_(self, *args,**kwds):
1922        """
1923        Evals legendre_P
1924        numerically with mpmath.
1925        EXAMPLES::
1926            sage: legendre_P(10,3).n(75)
1927            8.097453000000000000000e6
1928        """
1929        try:
1930            step_parent = kwds['parent']
1931        except KeyError:
1932            step_parent = parent(args[-1])
1933
1934        import sage.libs.mpmath.all as mpmath
1935
1936        try:
1937            precision = step_parent.prec()
1938        except AttributeError:
1939            precision = mpmath.mp.prec
1940
1941        return mpmath.call(mpmath.legendre,args[0],args[-1],prec = precision)
1942   
1943
1944    def _eval_special_values_(self,*args):
1945        """
1946        Special values known.
1947        EXAMPLES:
1948           
1949            sage: var('n')
1950            n
1951            sage: legendre_P(n,1)
1952            1
1953            sage: legendre_P(n,-1)
1954            (-1)^n
1955        """
1956        if args[-1] == 1:
1957            return 1
1958       
1959        if args[-1] == -1:
1960            return (-1)**args[0]
1961       
1962        if (args[-1] == 0):
1963            try:
1964                return (1+(-1)**args[0])/2*binomial(args[0],args[0]/2)/4**(args[0]/2)
1965            except TypeError:
1966                pass
1967
1968        raise ValueError("Value not found")
1969
1970    def _eval_numpy_(self, *args):
1971        #TODO: numpy_eval with help of the new scipy version!!!!
1972        #Reason scipy 0.8 supports stable and fast numerical evaluation
1973        #of ortho polys
1974        #Now this only evaluates the array pointwise, and only the first one...
1975        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
1976        """
1977        EXAMPLES::
1978        sage: import numpy
1979        sage: z = numpy.array([1,2])
1980        sage: z2 = numpy.array([[1,2],[1,2]])
1981        sage: z3 = numpy.array([1,2,3.])
1982        sage: legendre_P(1,z)
1983        array([1, 2])
1984        sage: legendre_P(1,z2)
1985        array([[1, 2],
1986              [1, 2]])
1987        sage: legendre_P(1,z3)
1988        array([ 1.,  2.,  3.])
1989
1990        """
1991        if isinstance(args[0],numpy.ndarray):
1992            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
1993
1994        result = numpy.zeros(args[-1].shape).tolist()
1995        if isinstance(args[-1][0],int):
1996            for k in xrange(len(args[-1])):
1997                result[k] = legendre_P(args[0],ZZ(args[-1][k]))
1998
1999        if isinstance(args[-1][0],float):
2000            for k in xrange(len(args[-1])):
2001                result[k] = legendre_P(args[0],RR(args[-1][k]))
2002
2003        if isinstance(args[-1][0],numpy.ndarray):
2004            for k in xrange(len(args[-1])):
2005                result[k] = legendre_P(args[0],args[-1][k])
2006
2007        return numpy.array(result)
2008
2009    def _derivative_(self,*args,**kwds):
2010        """return the derivative of legendre_P in
2011        form of the Legendre Polynomial.
2012        EXAMPLES::
2013        sage: n = var('n')
2014        sage: derivative(legendre_P(n,x),x)
2015        (n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x^2 - 1)
2016        sage: derivative(legendre_P(3,x),x)
2017        15/2*x^2 - 3/2
2018        """
2019        diff_param = kwds['diff_param']
2020        if diff_param == 0: 
2021            raise NotImplementedError(
2022"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
2023        else:
2024            return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*legendre_P(args[0],args[-1]))/(1-args[-1]**2)
2025       
2026
2027legendre_P = Func_legendre_P()     
2028
2029class Func_legendre_Q(OrthogonalPolynomial):
2030    """
2031    Return the chebyshev function of the second kind.
2032   
2033    REFERENCE:
2034   
2035    - A.S. 8 (p. 332)
2036   
2037    EXAMPLES::
2038   
2039        sage: P.<t> = QQ[]
2040        sage: legendre_Q(2, t)
2041        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
2042        sage: legendre_Q(3, 0.5)
2043        -0.198654771479482
2044        sage: legendre_Q(4, 2)
2045        443/16*I*pi + 443/16*log(3) - 365/12
2046        sage: legendre_Q(4, 2.0)
2047        0.00116107583162041 + 86.9828465962674*I
2048
2049    """
2050    def __init__(self):
2051        OrthogonalPolynomial.__init__(self,"legendre_Q",nargs = 2,
2052conversions =dict(maxima='legendre_q',mathematica='LegendreQ'))
2053
2054    def _maxima_init_evaled_(self, *args):
2055        """
2056        Maxima seems just fine for legendre Q. So we use it here!
2057        """
2058        n = args[0]
2059        x = args[1]
2060        return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
2061
2062    def _clenshaw_method_(self,*args):
2063        """
2064        Clenshaw method for legendre_q (means use the recursion...)
2065        This is much faster for numerical evaluation than maxima!
2066        See A.S. 8.5.3 (p. 334) for details for the recurions.
2067        Warning: The clanshaw method for the Legendre fUNCTIONS
2068        should only used for exact data types, when high orders are
2069        used, due to weak instabilities of the recursion!
2070        """
2071        raise NotImplementedError("Function not ready yet...")
2072       
2073        k = args[0]
2074        x = args[-1]
2075       
2076        if k == 0:
2077            return ln((1+x)/(1-x))/2
2078        elif k == 1:
2079            return x/2*ln((1+x)/(1-x))-1
2080        else:
2081            if is_Expression(x):
2082                raise NotImplementedError("Maxima works fine here!") 
2083                #it seems that the old method just works fine here...
2084                #raise NotImplementedError("clenshaw does not work well...")
2085            else:
2086                help1 = ln((1+x)/(1-x))/2
2087                help2 = x/2*ln((1+x)/(1-x))-1
2088
2089                for j in xrange(1,k):
2090                    help3 = (2*j+1)*x*help2 - j*help1
2091                    help3 = help3/(j+1)
2092                    help1 = help2
2093                    help2 = help3
2094               
2095            return help3   
2096
2097    def _eval_special_values_(self,*args):
2098        """
2099        Special values known.
2100        EXAMPLES:
2101           
2102            sage: var('n')
2103            n
2104            sage: legendre_Q(n,0)
2105            -1/2*sqrt(pi)*gamma(1/2*n + 1/2)*sin(1/2*pi*n)/gamma(1/2*n + 1)
2106        """
2107        if args[-1] == 1:
2108            return NaN
2109       
2110        if args[-1] == -1:
2111            return NaN
2112       
2113        if (args[-1] == 0):
2114            if is_Expression(args[0]): 
2115                try:
2116                    return -(sqrt(SR.pi()))/2*sin(SR.pi()/2*args[0])*gamma((args[0]+1)/2)/gamma(args[0]/2 + 1)
2117                except TypeError:
2118                    pass
2119            else:
2120                return -(sqrt(math.pi))/2*sin(math.pi/2*args[0])*gamma((args[0]+1)/2)/gamma(args[0]/2. + 1)
2121
2122        raise ValueError("Value not found")
2123
2124    def _evalf_(self, *args,**kwds):
2125        """
2126        Evals legendre_Q
2127        numerically with mpmath.
2128        EXAMPLES::
2129            sage: legendre_Q(10,3).n(75)
2130            2.079454941572578263731e-9 + 1.271944942879431601408e7*I
2131        """
2132
2133        try:
2134            step_parent = kwds['parent']
2135        except KeyError:
2136            step_parent = parent(args[-1])
2137
2138        import sage.libs.mpmath.all as mpmath
2139
2140        try:
2141            precision = step_parent.prec()
2142        except AttributeError:
2143            precision = mpmath.mp.prec
2144
2145        return conjugate(mpmath.call(mpmath.legenq,args[0],0,args[-1],prec = precision))
2146        #it seems that mpmath uses here a different branch of the logarithm
2147
2148    def _eval_numpy_(self, *args):
2149        #TODO: numpy_eval with help of the new scipy version!!!!
2150        #Reason scipy 0.8 supports stable and fast numerical evaluation
2151        #of ortho polys
2152        #Now this only evaluates the array pointwise, and only the first one...
2153        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
2154        """
2155        EXAMPLES::
2156        sage: import numpy
2157        sage: z = numpy.array([1,2])
2158        sage: z2 = numpy.array([[1,2],[1,2]])
2159        sage: z3 = numpy.array([1,2,3.])
2160        sage: legendre_Q(1,z/5.)
2161        array([-0.95945349, -0.83054043])
2162        sage: legendre_Q(1,z2/5.)
2163        array([[-0.95945349, -0.83054043],
2164               [-0.95945349, -0.83054043]])
2165        sage: legendre_Q(1,z3/5.)
2166        array([-0.95945349, -0.83054043, -0.58411169])
2167
2168        """
2169        if isinstance(args[0],numpy.ndarray):
2170            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
2171
2172        result = numpy.zeros(args[-1].shape).tolist()
2173        if isinstance(args[-1][0],int):
2174            for k in xrange(len(args[-1])):
2175                result[k] = legendre_Q(args[0],ZZ(args[-1][k]))
2176
2177        if isinstance(args[-1][0],float):
2178            for k in xrange(len(args[-1])):
2179                result[k] = legendre_Q(args[0],RR(args[-1][k]))
2180
2181        if isinstance(args[-1][0],numpy.ndarray):
2182            for k in xrange(len(args[-1])):
2183                result[k] = legendre_Q(args[0],args[-1][k])
2184
2185        return numpy.array(result)
2186
2187    def _derivative_(self,*args,**kwds):
2188        """return the derivative of legendre_Q in
2189        form of the Legendre Function.
2190        EXAMPLES::
2191        n = var('n')
2192        derivative(legendre_Q(n,x),x)
2193        (n*x*legendre_Q(n, x) - n*legendre_Q(n - 1, x))/(x^2 - 1)
2194        sage: derivative(legendre_Q(3,x),x)
2195        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))
2196        """
2197        diff_param = kwds['diff_param']
2198        if diff_param == 0: 
2199            raise NotImplementedError(
2200"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
2201        else:
2202            return (args[0]*args[-1]*legendre_Q(args[0],args[-1])-args[0]*legendre_Q(args[0]-1,args[-1]))/(args[-1]**2-1)
2203
2204legendre_Q = Func_legendre_Q()
2205
2206
2207class Func_gen_legendre_P(OrthogonalPolynomial):
2208
2209    def __init__(self):
2210        OrthogonalPolynomial.__init__(self,"gen_legendre_P",nargs = 3,
2211conversions =dict(maxima='assoc_legendre_p',mathematica='LegendreP'))
2212
2213    def _evalf_(self, *args,**kwds):
2214        """
2215        Evals gen_legendre_P
2216        numerically with mpmath.
2217        EXAMPLES::
2218            sage: gen_legendre_P(10,2,3).n(75)
2219            -7.194963600000000000000e8
2220        """
2221
2222        try:
2223            step_parent = kwds['parent']
2224        except KeyError:
2225            step_parent = parent(args[-1])
2226
2227        import sage.libs.mpmath.all as mpmath
2228
2229        try:
2230            precision = step_parent.prec()
2231        except AttributeError:
2232            precision = mpmath.mp.prec
2233
2234        return mpmath.call(mpmath.legenp,args[0],args[1],args[-1],prec = precision)
2235
2236    def _eval_special_values_(self,*args):
2237        """
2238        Special values known.
2239        EXAMPLES:
2240           
2241            sage: n, m = var('n m')
2242            sage: gen_legendre_P(n,m,0)
2243            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))
2244        """
2245           
2246        if args[1] == 0:
2247            return legendre_P(args[0],args[-1])
2248
2249        if (args[-1] == 0):
2250            if is_Expression(args[0]): 
2251                try:
2252                    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])
2253                except TypeError:
2254                    pass
2255            else:
2256                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]
2257
2258        raise ValueError("Value not found")
2259
2260    def _maxima_init_evaled_(self, *args):
2261        n = args[0]
2262        m = args[1]
2263        x = args[2]
2264        if is_Expression(n) or is_Expression(m):
2265            return None
2266       
2267        from sage.functions.all import sqrt
2268   
2269        if m.mod(2).is_zero() or m.is_one():
2270            return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
2271        else:
2272            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))
2273
2274
2275    def _eval_numpy_(self, *args):
2276        #TODO: numpy_eval with help of the new scipy version!!!!
2277        #Reason scipy 0.8 supports stable and fast numerical evaluation
2278        #of ortho polys
2279        #Now this only evaluates the array pointwise, and only the first one...
2280        #if isinstance(arg[0], numpy.ndarray). 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_P(1,1,z/5.)
2288        array([-0.9797959 , -0.91651514])
2289        sage: gen_legendre_P(1,1,z2/5.)
2290        array([[-0.9797959 , -0.91651514],
2291               [-0.9797959 , -0.91651514]])
2292        sage: gen_legendre_P(1,1,z3/5.)
2293        array([-0.9797959 , -0.91651514, -0.8       ])
2294        """
2295
2296        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
2297            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
2298       
2299        result = numpy.zeros(args[-1].shape).tolist()
2300        if isinstance(args[-1][0],int):
2301            for k in xrange(len(args[-1])):
2302                result[k] = gen_legendre_P(args[0],args[1],ZZ(args[-1][k]))
2303
2304        if isinstance(args[-1][0],float):
2305            for k in xrange(len(args[-1])):
2306                result[k] = gen_legendre_P(args[0],args[1],RR(args[-1][k]))
2307
2308        if isinstance(args[-1][0],numpy.ndarray):
2309            for k in xrange(len(args[-1])):
2310                result[k] = gen_legendre_P(args[0],args[1],args[-1][k])
2311
2312        return numpy.array(result)
2313     
2314    def _derivative_(self,*args,**kwds):
2315        """return the derivative of gen_legendre_P in
2316        form of the Legendre Function.
2317        EXAMPLES::
2318        sage: n,m = var('n m')
2319        sage: derivative(gen_legendre_P(n,m,x),x)
2320        (n*x*gen_legendre_P(n, m, x) - (m + n)*gen_legendre_P(n - 1, m, x))/(x^2 - 1)
2321        sage: derivative(gen_legendre_P(3,1,x),x)
2322        -15*sqrt(-x^2 + 1)*x + 3/2*(5*(x - 1)^2 + 10*x - 6)*x/sqrt(-x^2 + 1)
2323        """
2324        diff_param = kwds['diff_param']
2325        if diff_param in [0,1]: 
2326            raise NotImplementedError("Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
2327        else:
2328            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)
2329
2330gen_legendre_P = Func_gen_legendre_P()
2331
2332class Func_gen_legendre_Q(OrthogonalPolynomial):
2333
2334    def __init__(self):
2335        OrthogonalPolynomial.__init__(self,"gen_legendre_Q",nargs = 3,
2336conversions =dict(maxima='assoc_legendre_q',mathematica='LegendreQ'))
2337
2338    def _evalf_(self, *args,**kwds):
2339        """
2340        Evals gen_legendre_Q
2341        numerically with mpmath.
2342        EXAMPLES::
2343            sage: gen_legendre_Q(10,2,3).n(75)
2344            -2.773909528741569374688e-7 - 1.130182239430298584113e9*I
2345        """
2346
2347        try:
2348            step_parent = kwds['parent']
2349        except KeyError:
2350            step_parent = parent(args[-1])
2351
2352        import sage.libs.mpmath.all as mpmath
2353
2354        try:
2355            precision = step_parent.prec()
2356        except AttributeError:
2357            precision = mpmath.mp.prec
2358
2359        return mpmath.call(mpmath.legenq,args[0],args[1],args[-1],prec = precision)
2360
2361    def _eval_special_values_(self,*args):
2362        """
2363        Special values known.
2364        EXAMPLES:
2365           
2366            sage: n, m = var('n m')
2367            sage: gen_legendre_Q(n,m,0)
2368            -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)
2369        """
2370           
2371        if args[1] == 0:
2372            return legendre_Q(args[0],args[-1])
2373
2374        if (args[-1] == 0):
2375            if is_Expression(args[0]): 
2376                try:
2377                    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)
2378                except TypeError:
2379                    pass
2380            else:
2381                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]
2382
2383        raise ValueError("Value not found")
2384
2385    def _maxima_init_evaled_(self, *args):
2386        n = args[0]
2387        m = args[1]
2388        x = args[2]
2389        if is_Expression(n) or is_Expression(m):
2390            return None
2391       
2392        from sage.functions.all import sqrt
2393   
2394        if m <= n:
2395            return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
2396        if m == n + 1 or n == 0:
2397            if m.mod(2).is_zero():
2398                denom = (1 - x**2)**(m/2)
2399            else:
2400                denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)   
2401                if m == n + 1:
2402                    return (-1)**m*(m-1).factorial()*2**n/denom
2403                else:
2404                    return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
2405        else:
2406            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)
2407
2408
2409    def _eval_numpy_(self, *args):
2410        #TODO: numpy_eval with help of the new scipy version!!!!
2411        #Reason scipy 0.8 supports stable and fast numerical evaluation
2412        #of ortho polys
2413        #Now this only evaluates the array pointwise, and only the first one...
2414        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
2415        """
2416        EXAMPLES::
2417        sage: import numpy
2418        sage: z = numpy.array([1,2])
2419        sage: z2 = numpy.array([[1,2],[1,2]])
2420        sage: z3 = numpy.array([1,2,3.])
2421        sage: gen_legendre_Q(1,1,z/5.)
2422        array([-0.40276067, -0.82471644])
2423        sage: gen_legendre_Q(1,1,z2/5.)
2424        array([[-0.40276067, -0.82471644],
2425               [-0.40276067, -0.82471644]])
2426        sage: gen_legendre_Q(1,1,z3/5.)
2427        array([-0.40276067, -0.82471644, -1.30451774])
2428        """
2429
2430        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
2431            raise NotImplementedError("Support for numpy array in first argument(s) is not supported yet!")
2432       
2433        result = numpy.zeros(args[-1].shape).tolist()
2434        if isinstance(args[-1][0],int):
2435            for k in xrange(len(args[-1])):
2436                result[k] = gen_legendre_Q(args[0],args[1],ZZ(args[-1][k]))
2437
2438        if isinstance(args[-1][0],float):
2439            for k in xrange(len(args[-1])):
2440                result[k] = gen_legendre_Q(args[0],args[1],RR(args[-1][k]))
2441
2442        if isinstance(args[-1][0],numpy.ndarray):
2443            for k in xrange(len(args[-1])):
2444                result[k] = gen_legendre_Q(args[0],args[1],args[-1][k])
2445
2446        return numpy.array(result)
2447     
2448    def _derivative_(self,*args,**kwds):
2449        """return the derivative of gen_legendre_Q in
2450        form of the Legendre Function.
2451        EXAMPLES::
2452        sage: n,m = var('n m')
2453        sage: derivative(gen_legendre_Q(n,m,x),x)
2454        (n*x*gen_legendre_Q(n, m, x) - (m + n)*gen_legendre_Q(n - 1, m, x))/(x^2 - 1)
2455        sage: derivative(gen_legendre_Q(0,1,x),x)
2456        -x/(-x^2 + 1)^(3/2)
2457        """
2458        diff_param = kwds['diff_param']
2459        if diff_param in [0,1]: 
2460            raise NotImplementedError("Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
2461        else:
2462            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)
2463
2464
2465gen_legendre_Q = Func_gen_legendre_Q()
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476
2477