source: sage/functions/transcendental.py @ 3291:69aa0e79006c

Revision 3291:69aa0e79006c, 6.5 KB checked in by Carl Witty <cwitty@…>, 6 years ago (diff)

Fix one last doctest (for 64-bit)
Line 
1"""
2Transcendental Functions
3"""
4
5#*****************************************************************************
6#       Copyright (C) 2005 William Stein <wstein@gmail.com>
7#
8#  Distributed under the terms of the GNU General Public License (GPL)
9#
10#    This code is distributed in the hope that it will be useful,
11#    but WITHOUT ANY WARRANTY; without even the implied warranty of
12#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13#    General Public License for more details.
14#
15#  The full text of the GPL is available at:
16#
17#                  http://www.gnu.org/licenses/
18#*****************************************************************************
19
20import  sage.libs.pari.all
21pari = sage.libs.pari.all.pari
22import sage.rings.complex_field as complex_field
23import sage.rings.real_mpfr as real_field
24import sage.rings.complex_number
25
26from sage.rings.all import is_RealNumber, RealField, is_ComplexNumber, ComplexField
27
28def __prep_num(x):
29    if isinstance(x, sage.rings.complex_number.ComplexNumber):
30        x = str(x).replace("i","I")
31    else:
32        x = str(x)
33    return x
34
35CC = complex_field.ComplexField()
36I = CC.gen(0)
37def __eval(x):
38    return eval(x)
39
40def exponential_integral_1(x, n=0):
41    r"""
42    Returns the exponential integral $E_1(x)$. If the optional argument
43    $n$ is given, computes list of the first $n$ values of the exponential
44    integral $E_1(x m)$.
45
46    The exponential integral $E_1(x)$ is
47    $$
48             E_1(x) = \int_{x}^{\infty} e^{-t}/t dt
49    $$
50   
51    INPUT:
52        x -- a positive real number
53
54        n -- (default: 0) a nonnegative integer; if nonzero,
55             then return a list of values E_1(x*m) for
56             m = 1,2,3,...,n.   This is useful, e.g., when
57             computing derivatives of L-functions.
58
59    OUTPUT:
60        float -- if n is 0 (the default)
61      or
62        list -- list of floats if n > 0 
63
64    EXAMPLES:
65        sage: exponential_integral_1(2)
66        0.048900510708061118
67        sage: w = exponential_integral_1(2,4); w
68        [0.048900510708061118, 0.0037793524098489067, 0.00036008245216265867, 3.76656228439249e-05]
69
70
71    IMPLEMENTATION: We use the PARI C-library functions eint1 and
72    veceint1. 
73       
74    REFERENCE: See page 262, Prop 5.6.12, of Cohen's book "A Course
75    in Computational Algebraic Number Theory".
76
77    REMARKS: When called with the optional argument n, the PARI
78    C-library is fast for values of n up to some bound, then very
79    very slow.  For example, if x=5, then the computation takes less
80    than a second for n=800000, and takes "forever" for n=900000.
81   
82    """
83    if n <= 0:
84        return float(pari(x).eint1())
85    else:
86        return [float(z) for z in pari(x).eint1(n)]
87
88def gamma(s):
89    """
90    Gamma function at s.
91
92    EXAMPLES:
93        sage: gamma(CDF(0.5,14))
94        -4.05370307804e-10 - 5.77329961615e-10*I
95        sage: gamma(I)
96        -0.154949828301811 - 0.498015668118356*I
97        sage: gamma(6)
98        120.000000000000
99    """
100    try:
101        return s.gamma()
102    except AttributeError:
103        return CC(s).gamma()
104
105def gamma_inc(s, t):
106    """
107    Incomplete Gamma function Gamma(s,t).
108
109    EXAMPLES:
110        sage: gamma_inc(CDF(0,1), 3)
111        0.00320857499337 + 0.0124061862007*I
112        sage: gamma_inc(3, 3)
113        0.846380162253687 + 2.52435489670724e-29*I      # 32-bit
114        0.846380162253687                               # 64-bit
115        sage: gamma_inc(RDF(1), 3)
116        0.0497870683678639
117    """
118    try:
119        return s.gamma_inc(t)
120    except AttributeError:
121        if not (is_ComplexNumber(s)):
122            if is_ComplexNumber(t):
123                C = t.parent()
124            else:
125                C = ComplexField()
126            s = C(s)
127        return s.gamma_inc(t)
128
129
130# synonym.
131incomplete_gamma = gamma_inc
132
133def zeta(s):
134    """
135    Riemann zeta function at s with s a real or complex number.
136
137    INPUT:
138        s -- real or complex number
139
140    If s is a real number the computation is done using the MPFR
141    library.  When the input is not real, the computation is done
142    using the PARI C library.
143   
144    EXAMPLES:
145        sage: zeta(2)
146        1.64493406684823
147        sage: RR = RealField(200)
148        sage: zeta(RR(2))
149        1.6449340668482264364724151666460251892189499012067984377356
150    """
151    try:
152        return s.zeta()
153    except AttributeError:
154        return RealField()(s).zeta()
155
156##     prec = s.prec()
157##     s = pari.new_with_prec(s, prec)
158##     z = s.zeta()._sage_()
159##     if z.prec() < prec:
160##         raise RuntimeError, "Error computing zeta(%s) -- precision loss."%s
161##     return z
162
163def zeta_symmetric(s):
164    r"""
165    Completed function $\xi(s)$ that satisfies $\xi(s) = \xi(1-s)$ and
166    has zeros at the same points as the Riemann zeta function.
167
168    INPUT:
169        s -- real or complex number
170
171    If s is a real number the computation is done using the MPFR
172    library.  When the input is not real, the computation is done
173    using the PARI C library.
174   
175    More precisely,
176    $$
177       xi(s) = \gamma(s/2 + 1) * (s-1) * \pi^{-s/2} * \zeta(s).
178    $$
179   
180    EXAMPLES:
181        sage: zeta_symmetric(0.7)
182        0.497580414651127
183        sage: zeta_symmetric(1-0.7)
184        0.497580414651127
185        sage: RR = RealField(200)
186        sage: zeta_symmetric(RR('0.7'))
187        0.49758041465112690357779107525638385212657443284080589766062
188        sage: I = CC.0
189        sage: zeta_symmetric(RR('0.5') + I*RR('14.0'))
190        0.000201294444235258 + 4.74338450462408e-20*I
191        sage: zeta_symmetric(RR('0.5') + I*RR('14.1'))
192        0.0000489893483255687 + 1.18584612615602e-20*I
193        sage: zeta_symmetric(RR('0.5') + I*RR('14.2'))
194        -0.0000868931282620101 - 2.03287907341032e-20*I
195
196    REFERENCE:
197      I copied the definition of xi from
198        \url{http://www.math.ubc.ca/~pugh/RiemannZeta/RiemannZetaLong.html}
199    """
200    if not (is_ComplexNumber(s) or is_RealNumber(s)):
201        s = RealField()(s)
202       
203    if s == 1:  # deal with poles, hopefully
204        return s.parent()(1/2)
205   
206    R = s.parent()
207    return (s/2 + 1).gamma()   *    (s-1)   * (R.pi()**(-s/2))  *  s.zeta()
208   
209
210##     # Use PARI on complex nubmer
211##     prec = s.prec()
212##     s = pari.new_with_bits_prec(s, prec)
213##     pi = pari.pi()
214##     w = (s/2 + 1).gamma() * (s-1) * pi **(-s/2) * s.zeta()
215##     z = w._sage_()
216##     if z.prec() < prec:
217##         raise RuntimeError, "Error computing zeta_symmetric(%s) -- precision loss."%s
218##     return z
219
220
221#def pi_approx(prec=53):
222#    """
223#    Return pi computed to prec bits of precision.
224#    """
225#   return real_field.RealField(prec).pi()
Note: See TracBrowser for help on using the repository browser.