Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/gamma_functions.py: 20%
579 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from math import prod
3from sympy.core import Add, S, Dummy, expand_func
4from sympy.core.expr import Expr
5from sympy.core.function import Function, ArgumentIndexError, PoleError
6from sympy.core.logic import fuzzy_and, fuzzy_not
7from sympy.core.numbers import Rational, pi, oo, I
8from sympy.core.power import Pow
9from sympy.functions.special.zeta_functions import zeta
10from sympy.functions.special.error_functions import erf, erfc, Ei
11from sympy.functions.elementary.complexes import re, unpolarify
12from sympy.functions.elementary.exponential import exp, log
13from sympy.functions.elementary.integers import ceiling, floor
14from sympy.functions.elementary.miscellaneous import sqrt
15from sympy.functions.elementary.trigonometric import sin, cos, cot
16from sympy.functions.combinatorial.numbers import bernoulli, harmonic
17from sympy.functions.combinatorial.factorials import factorial, rf, RisingFactorial
18from sympy.utilities.misc import as_int
20from mpmath import mp, workprec
21from mpmath.libmp.libmpf import prec_to_dps
23def intlike(n):
24 try:
25 as_int(n, strict=False)
26 return True
27 except ValueError:
28 return False
30###############################################################################
31############################ COMPLETE GAMMA FUNCTION ##########################
32###############################################################################
34class gamma(Function):
35 r"""
36 The gamma function
38 .. math::
39 \Gamma(x) := \int^{\infty}_{0} t^{x-1} e^{-t} \mathrm{d}t.
41 Explanation
42 ===========
44 The ``gamma`` function implements the function which passes through the
45 values of the factorial function (i.e., $\Gamma(n) = (n - 1)!$ when n is
46 an integer). More generally, $\Gamma(z)$ is defined in the whole complex
47 plane except at the negative integers where there are simple poles.
49 Examples
50 ========
52 >>> from sympy import S, I, pi, gamma
53 >>> from sympy.abc import x
55 Several special values are known:
57 >>> gamma(1)
58 1
59 >>> gamma(4)
60 6
61 >>> gamma(S(3)/2)
62 sqrt(pi)/2
64 The ``gamma`` function obeys the mirror symmetry:
66 >>> from sympy import conjugate
67 >>> conjugate(gamma(x))
68 gamma(conjugate(x))
70 Differentiation with respect to $x$ is supported:
72 >>> from sympy import diff
73 >>> diff(gamma(x), x)
74 gamma(x)*polygamma(0, x)
76 Series expansion is also supported:
78 >>> from sympy import series
79 >>> series(gamma(x), x, 0, 3)
80 1/x - EulerGamma + x*(EulerGamma**2/2 + pi**2/12) + x**2*(-EulerGamma*pi**2/12 - zeta(3)/3 - EulerGamma**3/6) + O(x**3)
82 We can numerically evaluate the ``gamma`` function to arbitrary precision
83 on the whole complex plane:
85 >>> gamma(pi).evalf(40)
86 2.288037795340032417959588909060233922890
87 >>> gamma(1+I).evalf(20)
88 0.49801566811835604271 - 0.15494982830181068512*I
90 See Also
91 ========
93 lowergamma: Lower incomplete gamma function.
94 uppergamma: Upper incomplete gamma function.
95 polygamma: Polygamma function.
96 loggamma: Log Gamma function.
97 digamma: Digamma function.
98 trigamma: Trigamma function.
99 beta: Euler Beta function.
101 References
102 ==========
104 .. [1] https://en.wikipedia.org/wiki/Gamma_function
105 .. [2] https://dlmf.nist.gov/5
106 .. [3] https://mathworld.wolfram.com/GammaFunction.html
107 .. [4] https://functions.wolfram.com/GammaBetaErf/Gamma/
109 """
111 unbranched = True
112 _singularities = (S.ComplexInfinity,)
114 def fdiff(self, argindex=1):
115 if argindex == 1:
116 return self.func(self.args[0])*polygamma(0, self.args[0])
117 else:
118 raise ArgumentIndexError(self, argindex)
120 @classmethod
121 def eval(cls, arg):
122 if arg.is_Number:
123 if arg is S.NaN:
124 return S.NaN
125 elif arg is oo:
126 return oo
127 elif intlike(arg):
128 if arg.is_positive:
129 return factorial(arg - 1)
130 else:
131 return S.ComplexInfinity
132 elif arg.is_Rational:
133 if arg.q == 2:
134 n = abs(arg.p) // arg.q
136 if arg.is_positive:
137 k, coeff = n, S.One
138 else:
139 n = k = n + 1
141 if n & 1 == 0:
142 coeff = S.One
143 else:
144 coeff = S.NegativeOne
146 coeff *= prod(range(3, 2*k, 2))
148 if arg.is_positive:
149 return coeff*sqrt(pi) / 2**n
150 else:
151 return 2**n*sqrt(pi) / coeff
153 def _eval_expand_func(self, **hints):
154 arg = self.args[0]
155 if arg.is_Rational:
156 if abs(arg.p) > arg.q:
157 x = Dummy('x')
158 n = arg.p // arg.q
159 p = arg.p - n*arg.q
160 return self.func(x + n)._eval_expand_func().subs(x, Rational(p, arg.q))
162 if arg.is_Add:
163 coeff, tail = arg.as_coeff_add()
164 if coeff and coeff.q != 1:
165 intpart = floor(coeff)
166 tail = (coeff - intpart,) + tail
167 coeff = intpart
168 tail = arg._new_rawargs(*tail, reeval=False)
169 return self.func(tail)*RisingFactorial(tail, coeff)
171 return self.func(*self.args)
173 def _eval_conjugate(self):
174 return self.func(self.args[0].conjugate())
176 def _eval_is_real(self):
177 x = self.args[0]
178 if x.is_nonpositive and x.is_integer:
179 return False
180 if intlike(x) and x <= 0:
181 return False
182 if x.is_positive or x.is_noninteger:
183 return True
185 def _eval_is_positive(self):
186 x = self.args[0]
187 if x.is_positive:
188 return True
189 elif x.is_noninteger:
190 return floor(x).is_even
192 def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):
193 return exp(loggamma(z))
195 def _eval_rewrite_as_factorial(self, z, **kwargs):
196 return factorial(z - 1)
198 def _eval_nseries(self, x, n, logx, cdir=0):
199 x0 = self.args[0].limit(x, 0)
200 if not (x0.is_Integer and x0 <= 0):
201 return super()._eval_nseries(x, n, logx)
202 t = self.args[0] - x0
203 return (self.func(t + 1)/rf(self.args[0], -x0 + 1))._eval_nseries(x, n, logx)
205 def _eval_as_leading_term(self, x, logx=None, cdir=0):
206 arg = self.args[0]
207 x0 = arg.subs(x, 0)
209 if x0.is_integer and x0.is_nonpositive:
210 n = -x0
211 res = S.NegativeOne**n/self.func(n + 1)
212 return res/(arg + n).as_leading_term(x)
213 elif not x0.is_infinite:
214 return self.func(x0)
215 raise PoleError()
218###############################################################################
219################## LOWER and UPPER INCOMPLETE GAMMA FUNCTIONS #################
220###############################################################################
222class lowergamma(Function):
223 r"""
224 The lower incomplete gamma function.
226 Explanation
227 ===========
229 It can be defined as the meromorphic continuation of
231 .. math::
232 \gamma(s, x) := \int_0^x t^{s-1} e^{-t} \mathrm{d}t = \Gamma(s) - \Gamma(s, x).
234 This can be shown to be the same as
236 .. math::
237 \gamma(s, x) = \frac{x^s}{s} {}_1F_1\left({s \atop s+1} \middle| -x\right),
239 where ${}_1F_1$ is the (confluent) hypergeometric function.
241 Examples
242 ========
244 >>> from sympy import lowergamma, S
245 >>> from sympy.abc import s, x
246 >>> lowergamma(s, x)
247 lowergamma(s, x)
248 >>> lowergamma(3, x)
249 -2*(x**2/2 + x + 1)*exp(-x) + 2
250 >>> lowergamma(-S(1)/2, x)
251 -2*sqrt(pi)*erf(sqrt(x)) - 2*exp(-x)/sqrt(x)
253 See Also
254 ========
256 gamma: Gamma function.
257 uppergamma: Upper incomplete gamma function.
258 polygamma: Polygamma function.
259 loggamma: Log Gamma function.
260 digamma: Digamma function.
261 trigamma: Trigamma function.
262 beta: Euler Beta function.
264 References
265 ==========
267 .. [1] https://en.wikipedia.org/wiki/Incomplete_gamma_function#Lower_incomplete_gamma_function
268 .. [2] Abramowitz, Milton; Stegun, Irene A., eds. (1965), Chapter 6,
269 Section 5, Handbook of Mathematical Functions with Formulas, Graphs,
270 and Mathematical Tables
271 .. [3] https://dlmf.nist.gov/8
272 .. [4] https://functions.wolfram.com/GammaBetaErf/Gamma2/
273 .. [5] https://functions.wolfram.com/GammaBetaErf/Gamma3/
275 """
278 def fdiff(self, argindex=2):
279 from sympy.functions.special.hyper import meijerg
280 if argindex == 2:
281 a, z = self.args
282 return exp(-unpolarify(z))*z**(a - 1)
283 elif argindex == 1:
284 a, z = self.args
285 return gamma(a)*digamma(a) - log(z)*uppergamma(a, z) \
286 - meijerg([], [1, 1], [0, 0, a], [], z)
288 else:
289 raise ArgumentIndexError(self, argindex)
291 @classmethod
292 def eval(cls, a, x):
293 # For lack of a better place, we use this one to extract branching
294 # information. The following can be
295 # found in the literature (c/f references given above), albeit scattered:
296 # 1) For fixed x != 0, lowergamma(s, x) is an entire function of s
297 # 2) For fixed positive integers s, lowergamma(s, x) is an entire
298 # function of x.
299 # 3) For fixed non-positive integers s,
300 # lowergamma(s, exp(I*2*pi*n)*x) =
301 # 2*pi*I*n*(-1)**(-s)/factorial(-s) + lowergamma(s, x)
302 # (this follows from lowergamma(s, x).diff(x) = x**(s-1)*exp(-x)).
303 # 4) For fixed non-integral s,
304 # lowergamma(s, x) = x**s*gamma(s)*lowergamma_unbranched(s, x),
305 # where lowergamma_unbranched(s, x) is an entire function (in fact
306 # of both s and x), i.e.
307 # lowergamma(s, exp(2*I*pi*n)*x) = exp(2*pi*I*n*a)*lowergamma(a, x)
308 if x is S.Zero:
309 return S.Zero
310 nx, n = x.extract_branch_factor()
311 if a.is_integer and a.is_positive:
312 nx = unpolarify(x)
313 if nx != x:
314 return lowergamma(a, nx)
315 elif a.is_integer and a.is_nonpositive:
316 if n != 0:
317 return 2*pi*I*n*S.NegativeOne**(-a)/factorial(-a) + lowergamma(a, nx)
318 elif n != 0:
319 return exp(2*pi*I*n*a)*lowergamma(a, nx)
321 # Special values.
322 if a.is_Number:
323 if a is S.One:
324 return S.One - exp(-x)
325 elif a is S.Half:
326 return sqrt(pi)*erf(sqrt(x))
327 elif a.is_Integer or (2*a).is_Integer:
328 b = a - 1
329 if b.is_positive:
330 if a.is_integer:
331 return factorial(b) - exp(-x) * factorial(b) * Add(*[x ** k / factorial(k) for k in range(a)])
332 else:
333 return gamma(a)*(lowergamma(S.Half, x)/sqrt(pi) - exp(-x)*Add(*[x**(k - S.Half)/gamma(S.Half + k) for k in range(1, a + S.Half)]))
335 if not a.is_Integer:
336 return S.NegativeOne**(S.Half - a)*pi*erf(sqrt(x))/gamma(1 - a) + exp(-x)*Add(*[x**(k + a - 1)*gamma(a)/gamma(a + k) for k in range(1, Rational(3, 2) - a)])
338 if x.is_zero:
339 return S.Zero
341 def _eval_evalf(self, prec):
342 if all(x.is_number for x in self.args):
343 a = self.args[0]._to_mpmath(prec)
344 z = self.args[1]._to_mpmath(prec)
345 with workprec(prec):
346 res = mp.gammainc(a, 0, z)
347 return Expr._from_mpmath(res, prec)
348 else:
349 return self
351 def _eval_conjugate(self):
352 x = self.args[1]
353 if x not in (S.Zero, S.NegativeInfinity):
354 return self.func(self.args[0].conjugate(), x.conjugate())
356 def _eval_is_meromorphic(self, x, a):
357 # By https://en.wikipedia.org/wiki/Incomplete_gamma_function#Holomorphic_extension,
358 # lowergamma(s, z) = z**s*gamma(s)*gammastar(s, z),
359 # where gammastar(s, z) is holomorphic for all s and z.
360 # Hence the singularities of lowergamma are z = 0 (branch
361 # point) and nonpositive integer values of s (poles of gamma(s)).
362 s, z = self.args
363 args_merom = fuzzy_and([z._eval_is_meromorphic(x, a),
364 s._eval_is_meromorphic(x, a)])
365 if not args_merom:
366 return args_merom
367 z0 = z.subs(x, a)
368 if s.is_integer:
369 return fuzzy_and([s.is_positive, z0.is_finite])
370 s0 = s.subs(x, a)
371 return fuzzy_and([s0.is_finite, z0.is_finite, fuzzy_not(z0.is_zero)])
373 def _eval_aseries(self, n, args0, x, logx):
374 from sympy.series.order import O
375 s, z = self.args
376 if args0[0] is oo and not z.has(x):
377 coeff = z**s*exp(-z)
378 sum_expr = sum(z**k/rf(s, k + 1) for k in range(n - 1))
379 o = O(z**s*s**(-n))
380 return coeff*sum_expr + o
381 return super()._eval_aseries(n, args0, x, logx)
383 def _eval_rewrite_as_uppergamma(self, s, x, **kwargs):
384 return gamma(s) - uppergamma(s, x)
386 def _eval_rewrite_as_expint(self, s, x, **kwargs):
387 from sympy.functions.special.error_functions import expint
388 if s.is_integer and s.is_nonpositive:
389 return self
390 return self.rewrite(uppergamma).rewrite(expint)
392 def _eval_is_zero(self):
393 x = self.args[1]
394 if x.is_zero:
395 return True
398class uppergamma(Function):
399 r"""
400 The upper incomplete gamma function.
402 Explanation
403 ===========
405 It can be defined as the meromorphic continuation of
407 .. math::
408 \Gamma(s, x) := \int_x^\infty t^{s-1} e^{-t} \mathrm{d}t = \Gamma(s) - \gamma(s, x).
410 where $\gamma(s, x)$ is the lower incomplete gamma function,
411 :class:`lowergamma`. This can be shown to be the same as
413 .. math::
414 \Gamma(s, x) = \Gamma(s) - \frac{x^s}{s} {}_1F_1\left({s \atop s+1} \middle| -x\right),
416 where ${}_1F_1$ is the (confluent) hypergeometric function.
418 The upper incomplete gamma function is also essentially equivalent to the
419 generalized exponential integral:
421 .. math::
422 \operatorname{E}_{n}(x) = \int_{1}^{\infty}{\frac{e^{-xt}}{t^n} \, dt} = x^{n-1}\Gamma(1-n,x).
424 Examples
425 ========
427 >>> from sympy import uppergamma, S
428 >>> from sympy.abc import s, x
429 >>> uppergamma(s, x)
430 uppergamma(s, x)
431 >>> uppergamma(3, x)
432 2*(x**2/2 + x + 1)*exp(-x)
433 >>> uppergamma(-S(1)/2, x)
434 -2*sqrt(pi)*erfc(sqrt(x)) + 2*exp(-x)/sqrt(x)
435 >>> uppergamma(-2, x)
436 expint(3, x)/x**2
438 See Also
439 ========
441 gamma: Gamma function.
442 lowergamma: Lower incomplete gamma function.
443 polygamma: Polygamma function.
444 loggamma: Log Gamma function.
445 digamma: Digamma function.
446 trigamma: Trigamma function.
447 beta: Euler Beta function.
449 References
450 ==========
452 .. [1] https://en.wikipedia.org/wiki/Incomplete_gamma_function#Upper_incomplete_gamma_function
453 .. [2] Abramowitz, Milton; Stegun, Irene A., eds. (1965), Chapter 6,
454 Section 5, Handbook of Mathematical Functions with Formulas, Graphs,
455 and Mathematical Tables
456 .. [3] https://dlmf.nist.gov/8
457 .. [4] https://functions.wolfram.com/GammaBetaErf/Gamma2/
458 .. [5] https://functions.wolfram.com/GammaBetaErf/Gamma3/
459 .. [6] https://en.wikipedia.org/wiki/Exponential_integral#Relation_with_other_functions
461 """
464 def fdiff(self, argindex=2):
465 from sympy.functions.special.hyper import meijerg
466 if argindex == 2:
467 a, z = self.args
468 return -exp(-unpolarify(z))*z**(a - 1)
469 elif argindex == 1:
470 a, z = self.args
471 return uppergamma(a, z)*log(z) + meijerg([], [1, 1], [0, 0, a], [], z)
472 else:
473 raise ArgumentIndexError(self, argindex)
475 def _eval_evalf(self, prec):
476 if all(x.is_number for x in self.args):
477 a = self.args[0]._to_mpmath(prec)
478 z = self.args[1]._to_mpmath(prec)
479 with workprec(prec):
480 res = mp.gammainc(a, z, mp.inf)
481 return Expr._from_mpmath(res, prec)
482 return self
484 @classmethod
485 def eval(cls, a, z):
486 from sympy.functions.special.error_functions import expint
487 if z.is_Number:
488 if z is S.NaN:
489 return S.NaN
490 elif z is oo:
491 return S.Zero
492 elif z.is_zero:
493 if re(a).is_positive:
494 return gamma(a)
496 # We extract branching information here. C/f lowergamma.
497 nx, n = z.extract_branch_factor()
498 if a.is_integer and a.is_positive:
499 nx = unpolarify(z)
500 if z != nx:
501 return uppergamma(a, nx)
502 elif a.is_integer and a.is_nonpositive:
503 if n != 0:
504 return -2*pi*I*n*S.NegativeOne**(-a)/factorial(-a) + uppergamma(a, nx)
505 elif n != 0:
506 return gamma(a)*(1 - exp(2*pi*I*n*a)) + exp(2*pi*I*n*a)*uppergamma(a, nx)
508 # Special values.
509 if a.is_Number:
510 if a is S.Zero and z.is_positive:
511 return -Ei(-z)
512 elif a is S.One:
513 return exp(-z)
514 elif a is S.Half:
515 return sqrt(pi)*erfc(sqrt(z))
516 elif a.is_Integer or (2*a).is_Integer:
517 b = a - 1
518 if b.is_positive:
519 if a.is_integer:
520 return exp(-z) * factorial(b) * Add(*[z**k / factorial(k)
521 for k in range(a)])
522 else:
523 return (gamma(a) * erfc(sqrt(z)) +
524 S.NegativeOne**(a - S(3)/2) * exp(-z) * sqrt(z)
525 * Add(*[gamma(-S.Half - k) * (-z)**k / gamma(1-a)
526 for k in range(a - S.Half)]))
527 elif b.is_Integer:
528 return expint(-b, z)*unpolarify(z)**(b + 1)
530 if not a.is_Integer:
531 return (S.NegativeOne**(S.Half - a) * pi*erfc(sqrt(z))/gamma(1-a)
532 - z**a * exp(-z) * Add(*[z**k * gamma(a) / gamma(a+k+1)
533 for k in range(S.Half - a)]))
535 if a.is_zero and z.is_positive:
536 return -Ei(-z)
538 if z.is_zero and re(a).is_positive:
539 return gamma(a)
541 def _eval_conjugate(self):
542 z = self.args[1]
543 if z not in (S.Zero, S.NegativeInfinity):
544 return self.func(self.args[0].conjugate(), z.conjugate())
546 def _eval_is_meromorphic(self, x, a):
547 return lowergamma._eval_is_meromorphic(self, x, a)
549 def _eval_rewrite_as_lowergamma(self, s, x, **kwargs):
550 return gamma(s) - lowergamma(s, x)
552 def _eval_rewrite_as_tractable(self, s, x, **kwargs):
553 return exp(loggamma(s)) - lowergamma(s, x)
555 def _eval_rewrite_as_expint(self, s, x, **kwargs):
556 from sympy.functions.special.error_functions import expint
557 return expint(1 - s, x)*x**s
560###############################################################################
561###################### POLYGAMMA and LOGGAMMA FUNCTIONS #######################
562###############################################################################
564class polygamma(Function):
565 r"""
566 The function ``polygamma(n, z)`` returns ``log(gamma(z)).diff(n + 1)``.
568 Explanation
569 ===========
571 It is a meromorphic function on $\mathbb{C}$ and defined as the $(n+1)$-th
572 derivative of the logarithm of the gamma function:
574 .. math::
575 \psi^{(n)} (z) := \frac{\mathrm{d}^{n+1}}{\mathrm{d} z^{n+1}} \log\Gamma(z).
577 For `n` not a nonnegative integer the generalization by Espinosa and Moll [5]_
578 is used:
580 .. math:: \psi(s,z) = \frac{\zeta'(s+1, z) + (\gamma + \psi(-s)) \zeta(s+1, z)}
581 {\Gamma(-s)}
583 Examples
584 ========
586 Several special values are known:
588 >>> from sympy import S, polygamma
589 >>> polygamma(0, 1)
590 -EulerGamma
591 >>> polygamma(0, 1/S(2))
592 -2*log(2) - EulerGamma
593 >>> polygamma(0, 1/S(3))
594 -log(3) - sqrt(3)*pi/6 - EulerGamma - log(sqrt(3))
595 >>> polygamma(0, 1/S(4))
596 -pi/2 - log(4) - log(2) - EulerGamma
597 >>> polygamma(0, 2)
598 1 - EulerGamma
599 >>> polygamma(0, 23)
600 19093197/5173168 - EulerGamma
602 >>> from sympy import oo, I
603 >>> polygamma(0, oo)
604 oo
605 >>> polygamma(0, -oo)
606 oo
607 >>> polygamma(0, I*oo)
608 oo
609 >>> polygamma(0, -I*oo)
610 oo
612 Differentiation with respect to $x$ is supported:
614 >>> from sympy import Symbol, diff
615 >>> x = Symbol("x")
616 >>> diff(polygamma(0, x), x)
617 polygamma(1, x)
618 >>> diff(polygamma(0, x), x, 2)
619 polygamma(2, x)
620 >>> diff(polygamma(0, x), x, 3)
621 polygamma(3, x)
622 >>> diff(polygamma(1, x), x)
623 polygamma(2, x)
624 >>> diff(polygamma(1, x), x, 2)
625 polygamma(3, x)
626 >>> diff(polygamma(2, x), x)
627 polygamma(3, x)
628 >>> diff(polygamma(2, x), x, 2)
629 polygamma(4, x)
631 >>> n = Symbol("n")
632 >>> diff(polygamma(n, x), x)
633 polygamma(n + 1, x)
634 >>> diff(polygamma(n, x), x, 2)
635 polygamma(n + 2, x)
637 We can rewrite ``polygamma`` functions in terms of harmonic numbers:
639 >>> from sympy import harmonic
640 >>> polygamma(0, x).rewrite(harmonic)
641 harmonic(x - 1) - EulerGamma
642 >>> polygamma(2, x).rewrite(harmonic)
643 2*harmonic(x - 1, 3) - 2*zeta(3)
644 >>> ni = Symbol("n", integer=True)
645 >>> polygamma(ni, x).rewrite(harmonic)
646 (-1)**(n + 1)*(-harmonic(x - 1, n + 1) + zeta(n + 1))*factorial(n)
648 See Also
649 ========
651 gamma: Gamma function.
652 lowergamma: Lower incomplete gamma function.
653 uppergamma: Upper incomplete gamma function.
654 loggamma: Log Gamma function.
655 digamma: Digamma function.
656 trigamma: Trigamma function.
657 beta: Euler Beta function.
659 References
660 ==========
662 .. [1] https://en.wikipedia.org/wiki/Polygamma_function
663 .. [2] https://mathworld.wolfram.com/PolygammaFunction.html
664 .. [3] https://functions.wolfram.com/GammaBetaErf/PolyGamma/
665 .. [4] https://functions.wolfram.com/GammaBetaErf/PolyGamma2/
666 .. [5] O. Espinosa and V. Moll, "A generalized polygamma function",
667 *Integral Transforms and Special Functions* (2004), 101-115.
669 """
671 @classmethod
672 def eval(cls, n, z):
673 if n is S.NaN or z is S.NaN:
674 return S.NaN
675 elif z is oo:
676 return oo if n.is_zero else S.Zero
677 elif z.is_Integer and z.is_nonpositive:
678 return S.ComplexInfinity
679 elif n is S.NegativeOne:
680 return loggamma(z) - log(2*pi) / 2
681 elif n.is_zero:
682 if z is -oo or z.extract_multiplicatively(I) in (oo, -oo):
683 return oo
684 elif z.is_Integer:
685 return harmonic(z-1) - S.EulerGamma
686 elif z.is_Rational:
687 # TODO n == 1 also can do some rational z
688 p, q = z.as_numer_denom()
689 # only expand for small denominators to avoid creating long expressions
690 if q <= 6:
691 return expand_func(polygamma(S.Zero, z, evaluate=False))
692 elif n.is_integer and n.is_nonnegative:
693 nz = unpolarify(z)
694 if z != nz:
695 return polygamma(n, nz)
696 if z.is_Integer:
697 return S.NegativeOne**(n+1) * factorial(n) * zeta(n+1, z)
698 elif z is S.Half:
699 return S.NegativeOne**(n+1) * factorial(n) * (2**(n+1)-1) * zeta(n+1)
701 def _eval_is_real(self):
702 if self.args[0].is_positive and self.args[1].is_positive:
703 return True
705 def _eval_is_complex(self):
706 z = self.args[1]
707 is_negative_integer = fuzzy_and([z.is_negative, z.is_integer])
708 return fuzzy_and([z.is_complex, fuzzy_not(is_negative_integer)])
710 def _eval_is_positive(self):
711 n, z = self.args
712 if n.is_positive:
713 if n.is_odd and z.is_real:
714 return True
715 if n.is_even and z.is_positive:
716 return False
718 def _eval_is_negative(self):
719 n, z = self.args
720 if n.is_positive:
721 if n.is_even and z.is_positive:
722 return True
723 if n.is_odd and z.is_real:
724 return False
726 def _eval_expand_func(self, **hints):
727 n, z = self.args
729 if n.is_Integer and n.is_nonnegative:
730 if z.is_Add:
731 coeff = z.args[0]
732 if coeff.is_Integer:
733 e = -(n + 1)
734 if coeff > 0:
735 tail = Add(*[Pow(
736 z - i, e) for i in range(1, int(coeff) + 1)])
737 else:
738 tail = -Add(*[Pow(
739 z + i, e) for i in range(int(-coeff))])
740 return polygamma(n, z - coeff) + S.NegativeOne**n*factorial(n)*tail
742 elif z.is_Mul:
743 coeff, z = z.as_two_terms()
744 if coeff.is_Integer and coeff.is_positive:
745 tail = [polygamma(n, z + Rational(
746 i, coeff)) for i in range(int(coeff))]
747 if n == 0:
748 return Add(*tail)/coeff + log(coeff)
749 else:
750 return Add(*tail)/coeff**(n + 1)
751 z *= coeff
753 if n == 0 and z.is_Rational:
754 p, q = z.as_numer_denom()
756 # Reference:
757 # Values of the polygamma functions at rational arguments, J. Choi, 2007
758 part_1 = -S.EulerGamma - pi * cot(p * pi / q) / 2 - log(q) + Add(
759 *[cos(2 * k * pi * p / q) * log(2 * sin(k * pi / q)) for k in range(1, q)])
761 if z > 0:
762 n = floor(z)
763 z0 = z - n
764 return part_1 + Add(*[1 / (z0 + k) for k in range(n)])
765 elif z < 0:
766 n = floor(1 - z)
767 z0 = z + n
768 return part_1 - Add(*[1 / (z0 - 1 - k) for k in range(n)])
770 if n == -1:
771 return loggamma(z) - log(2*pi) / 2
772 if n.is_integer is False or n.is_nonnegative is False:
773 s = Dummy("s")
774 dzt = zeta(s, z).diff(s).subs(s, n+1)
775 return (dzt + (S.EulerGamma + digamma(-n)) * zeta(n+1, z)) / gamma(-n)
777 return polygamma(n, z)
779 def _eval_rewrite_as_zeta(self, n, z, **kwargs):
780 if n.is_integer and n.is_positive:
781 return S.NegativeOne**(n + 1)*factorial(n)*zeta(n + 1, z)
783 def _eval_rewrite_as_harmonic(self, n, z, **kwargs):
784 if n.is_integer:
785 if n.is_zero:
786 return harmonic(z - 1) - S.EulerGamma
787 else:
788 return S.NegativeOne**(n+1) * factorial(n) * (zeta(n+1) - harmonic(z-1, n+1))
790 def _eval_as_leading_term(self, x, logx=None, cdir=0):
791 from sympy.series.order import Order
792 n, z = [a.as_leading_term(x) for a in self.args]
793 o = Order(z, x)
794 if n == 0 and o.contains(1/x):
795 logx = log(x) if logx is None else logx
796 return o.getn() * logx
797 else:
798 return self.func(n, z)
800 def fdiff(self, argindex=2):
801 if argindex == 2:
802 n, z = self.args[:2]
803 return polygamma(n + 1, z)
804 else:
805 raise ArgumentIndexError(self, argindex)
807 def _eval_aseries(self, n, args0, x, logx):
808 from sympy.series.order import Order
809 if args0[1] != oo or not \
810 (self.args[0].is_Integer and self.args[0].is_nonnegative):
811 return super()._eval_aseries(n, args0, x, logx)
812 z = self.args[1]
813 N = self.args[0]
815 if N == 0:
816 # digamma function series
817 # Abramowitz & Stegun, p. 259, 6.3.18
818 r = log(z) - 1/(2*z)
819 o = None
820 if n < 2:
821 o = Order(1/z, x)
822 else:
823 m = ceiling((n + 1)//2)
824 l = [bernoulli(2*k) / (2*k*z**(2*k)) for k in range(1, m)]
825 r -= Add(*l)
826 o = Order(1/z**n, x)
827 return r._eval_nseries(x, n, logx) + o
828 else:
829 # proper polygamma function
830 # Abramowitz & Stegun, p. 260, 6.4.10
831 # We return terms to order higher than O(x**n) on purpose
832 # -- otherwise we would not be able to return any terms for
833 # quite a long time!
834 fac = gamma(N)
835 e0 = fac + N*fac/(2*z)
836 m = ceiling((n + 1)//2)
837 for k in range(1, m):
838 fac = fac*(2*k + N - 1)*(2*k + N - 2) / ((2*k)*(2*k - 1))
839 e0 += bernoulli(2*k)*fac/z**(2*k)
840 o = Order(1/z**(2*m), x)
841 if n == 0:
842 o = Order(1/z, x)
843 elif n == 1:
844 o = Order(1/z**2, x)
845 r = e0._eval_nseries(z, n, logx) + o
846 return (-1 * (-1/z)**N * r)._eval_nseries(x, n, logx)
848 def _eval_evalf(self, prec):
849 if not all(i.is_number for i in self.args):
850 return
851 s = self.args[0]._to_mpmath(prec+12)
852 z = self.args[1]._to_mpmath(prec+12)
853 if mp.isint(z) and z <= 0:
854 return S.ComplexInfinity
855 with workprec(prec+12):
856 if mp.isint(s) and s >= 0:
857 res = mp.polygamma(s, z)
858 else:
859 zt = mp.zeta(s+1, z)
860 dzt = mp.zeta(s+1, z, 1)
861 res = (dzt + (mp.euler + mp.digamma(-s)) * zt) * mp.rgamma(-s)
862 return Expr._from_mpmath(res, prec)
865class loggamma(Function):
866 r"""
867 The ``loggamma`` function implements the logarithm of the
868 gamma function (i.e., $\log\Gamma(x)$).
870 Examples
871 ========
873 Several special values are known. For numerical integral
874 arguments we have:
876 >>> from sympy import loggamma
877 >>> loggamma(-2)
878 oo
879 >>> loggamma(0)
880 oo
881 >>> loggamma(1)
882 0
883 >>> loggamma(2)
884 0
885 >>> loggamma(3)
886 log(2)
888 And for symbolic values:
890 >>> from sympy import Symbol
891 >>> n = Symbol("n", integer=True, positive=True)
892 >>> loggamma(n)
893 log(gamma(n))
894 >>> loggamma(-n)
895 oo
897 For half-integral values:
899 >>> from sympy import S
900 >>> loggamma(S(5)/2)
901 log(3*sqrt(pi)/4)
902 >>> loggamma(n/2)
903 log(2**(1 - n)*sqrt(pi)*gamma(n)/gamma(n/2 + 1/2))
905 And general rational arguments:
907 >>> from sympy import expand_func
908 >>> L = loggamma(S(16)/3)
909 >>> expand_func(L).doit()
910 -5*log(3) + loggamma(1/3) + log(4) + log(7) + log(10) + log(13)
911 >>> L = loggamma(S(19)/4)
912 >>> expand_func(L).doit()
913 -4*log(4) + loggamma(3/4) + log(3) + log(7) + log(11) + log(15)
914 >>> L = loggamma(S(23)/7)
915 >>> expand_func(L).doit()
916 -3*log(7) + log(2) + loggamma(2/7) + log(9) + log(16)
918 The ``loggamma`` function has the following limits towards infinity:
920 >>> from sympy import oo
921 >>> loggamma(oo)
922 oo
923 >>> loggamma(-oo)
924 zoo
926 The ``loggamma`` function obeys the mirror symmetry
927 if $x \in \mathbb{C} \setminus \{-\infty, 0\}$:
929 >>> from sympy.abc import x
930 >>> from sympy import conjugate
931 >>> conjugate(loggamma(x))
932 loggamma(conjugate(x))
934 Differentiation with respect to $x$ is supported:
936 >>> from sympy import diff
937 >>> diff(loggamma(x), x)
938 polygamma(0, x)
940 Series expansion is also supported:
942 >>> from sympy import series
943 >>> series(loggamma(x), x, 0, 4).cancel()
944 -log(x) - EulerGamma*x + pi**2*x**2/12 - x**3*zeta(3)/3 + O(x**4)
946 We can numerically evaluate the ``loggamma`` function
947 to arbitrary precision on the whole complex plane:
949 >>> from sympy import I
950 >>> loggamma(5).evalf(30)
951 3.17805383034794561964694160130
952 >>> loggamma(I).evalf(20)
953 -0.65092319930185633889 - 1.8724366472624298171*I
955 See Also
956 ========
958 gamma: Gamma function.
959 lowergamma: Lower incomplete gamma function.
960 uppergamma: Upper incomplete gamma function.
961 polygamma: Polygamma function.
962 digamma: Digamma function.
963 trigamma: Trigamma function.
964 beta: Euler Beta function.
966 References
967 ==========
969 .. [1] https://en.wikipedia.org/wiki/Gamma_function
970 .. [2] https://dlmf.nist.gov/5
971 .. [3] https://mathworld.wolfram.com/LogGammaFunction.html
972 .. [4] https://functions.wolfram.com/GammaBetaErf/LogGamma/
974 """
975 @classmethod
976 def eval(cls, z):
977 if z.is_integer:
978 if z.is_nonpositive:
979 return oo
980 elif z.is_positive:
981 return log(gamma(z))
982 elif z.is_rational:
983 p, q = z.as_numer_denom()
984 # Half-integral values:
985 if p.is_positive and q == 2:
986 return log(sqrt(pi) * 2**(1 - p) * gamma(p) / gamma((p + 1)*S.Half))
988 if z is oo:
989 return oo
990 elif abs(z) is oo:
991 return S.ComplexInfinity
992 if z is S.NaN:
993 return S.NaN
995 def _eval_expand_func(self, **hints):
996 from sympy.concrete.summations import Sum
997 z = self.args[0]
999 if z.is_Rational:
1000 p, q = z.as_numer_denom()
1001 # General rational arguments (u + p/q)
1002 # Split z as n + p/q with p < q
1003 n = p // q
1004 p = p - n*q
1005 if p.is_positive and q.is_positive and p < q:
1006 k = Dummy("k")
1007 if n.is_positive:
1008 return loggamma(p / q) - n*log(q) + Sum(log((k - 1)*q + p), (k, 1, n))
1009 elif n.is_negative:
1010 return loggamma(p / q) - n*log(q) + pi*I*n - Sum(log(k*q - p), (k, 1, -n))
1011 elif n.is_zero:
1012 return loggamma(p / q)
1014 return self
1016 def _eval_nseries(self, x, n, logx=None, cdir=0):
1017 x0 = self.args[0].limit(x, 0)
1018 if x0.is_zero:
1019 f = self._eval_rewrite_as_intractable(*self.args)
1020 return f._eval_nseries(x, n, logx)
1021 return super()._eval_nseries(x, n, logx)
1023 def _eval_aseries(self, n, args0, x, logx):
1024 from sympy.series.order import Order
1025 if args0[0] != oo:
1026 return super()._eval_aseries(n, args0, x, logx)
1027 z = self.args[0]
1028 r = log(z)*(z - S.Half) - z + log(2*pi)/2
1029 l = [bernoulli(2*k) / (2*k*(2*k - 1)*z**(2*k - 1)) for k in range(1, n)]
1030 o = None
1031 if n == 0:
1032 o = Order(1, x)
1033 else:
1034 o = Order(1/z**n, x)
1035 # It is very inefficient to first add the order and then do the nseries
1036 return (r + Add(*l))._eval_nseries(x, n, logx) + o
1038 def _eval_rewrite_as_intractable(self, z, **kwargs):
1039 return log(gamma(z))
1041 def _eval_is_real(self):
1042 z = self.args[0]
1043 if z.is_positive:
1044 return True
1045 elif z.is_nonpositive:
1046 return False
1048 def _eval_conjugate(self):
1049 z = self.args[0]
1050 if z not in (S.Zero, S.NegativeInfinity):
1051 return self.func(z.conjugate())
1053 def fdiff(self, argindex=1):
1054 if argindex == 1:
1055 return polygamma(0, self.args[0])
1056 else:
1057 raise ArgumentIndexError(self, argindex)
1060class digamma(Function):
1061 r"""
1062 The ``digamma`` function is the first derivative of the ``loggamma``
1063 function
1065 .. math::
1066 \psi(x) := \frac{\mathrm{d}}{\mathrm{d} z} \log\Gamma(z)
1067 = \frac{\Gamma'(z)}{\Gamma(z) }.
1069 In this case, ``digamma(z) = polygamma(0, z)``.
1071 Examples
1072 ========
1074 >>> from sympy import digamma
1075 >>> digamma(0)
1076 zoo
1077 >>> from sympy import Symbol
1078 >>> z = Symbol('z')
1079 >>> digamma(z)
1080 polygamma(0, z)
1082 To retain ``digamma`` as it is:
1084 >>> digamma(0, evaluate=False)
1085 digamma(0)
1086 >>> digamma(z, evaluate=False)
1087 digamma(z)
1089 See Also
1090 ========
1092 gamma: Gamma function.
1093 lowergamma: Lower incomplete gamma function.
1094 uppergamma: Upper incomplete gamma function.
1095 polygamma: Polygamma function.
1096 loggamma: Log Gamma function.
1097 trigamma: Trigamma function.
1098 beta: Euler Beta function.
1100 References
1101 ==========
1103 .. [1] https://en.wikipedia.org/wiki/Digamma_function
1104 .. [2] https://mathworld.wolfram.com/DigammaFunction.html
1105 .. [3] https://functions.wolfram.com/GammaBetaErf/PolyGamma2/
1107 """
1108 def _eval_evalf(self, prec):
1109 z = self.args[0]
1110 nprec = prec_to_dps(prec)
1111 return polygamma(0, z).evalf(n=nprec)
1113 def fdiff(self, argindex=1):
1114 z = self.args[0]
1115 return polygamma(0, z).fdiff()
1117 def _eval_is_real(self):
1118 z = self.args[0]
1119 return polygamma(0, z).is_real
1121 def _eval_is_positive(self):
1122 z = self.args[0]
1123 return polygamma(0, z).is_positive
1125 def _eval_is_negative(self):
1126 z = self.args[0]
1127 return polygamma(0, z).is_negative
1129 def _eval_aseries(self, n, args0, x, logx):
1130 as_polygamma = self.rewrite(polygamma)
1131 args0 = [S.Zero,] + args0
1132 return as_polygamma._eval_aseries(n, args0, x, logx)
1134 @classmethod
1135 def eval(cls, z):
1136 return polygamma(0, z)
1138 def _eval_expand_func(self, **hints):
1139 z = self.args[0]
1140 return polygamma(0, z).expand(func=True)
1142 def _eval_rewrite_as_harmonic(self, z, **kwargs):
1143 return harmonic(z - 1) - S.EulerGamma
1145 def _eval_rewrite_as_polygamma(self, z, **kwargs):
1146 return polygamma(0, z)
1148 def _eval_as_leading_term(self, x, logx=None, cdir=0):
1149 z = self.args[0]
1150 return polygamma(0, z).as_leading_term(x)
1154class trigamma(Function):
1155 r"""
1156 The ``trigamma`` function is the second derivative of the ``loggamma``
1157 function
1159 .. math::
1160 \psi^{(1)}(z) := \frac{\mathrm{d}^{2}}{\mathrm{d} z^{2}} \log\Gamma(z).
1162 In this case, ``trigamma(z) = polygamma(1, z)``.
1164 Examples
1165 ========
1167 >>> from sympy import trigamma
1168 >>> trigamma(0)
1169 zoo
1170 >>> from sympy import Symbol
1171 >>> z = Symbol('z')
1172 >>> trigamma(z)
1173 polygamma(1, z)
1175 To retain ``trigamma`` as it is:
1177 >>> trigamma(0, evaluate=False)
1178 trigamma(0)
1179 >>> trigamma(z, evaluate=False)
1180 trigamma(z)
1183 See Also
1184 ========
1186 gamma: Gamma function.
1187 lowergamma: Lower incomplete gamma function.
1188 uppergamma: Upper incomplete gamma function.
1189 polygamma: Polygamma function.
1190 loggamma: Log Gamma function.
1191 digamma: Digamma function.
1192 beta: Euler Beta function.
1194 References
1195 ==========
1197 .. [1] https://en.wikipedia.org/wiki/Trigamma_function
1198 .. [2] https://mathworld.wolfram.com/TrigammaFunction.html
1199 .. [3] https://functions.wolfram.com/GammaBetaErf/PolyGamma2/
1201 """
1202 def _eval_evalf(self, prec):
1203 z = self.args[0]
1204 nprec = prec_to_dps(prec)
1205 return polygamma(1, z).evalf(n=nprec)
1207 def fdiff(self, argindex=1):
1208 z = self.args[0]
1209 return polygamma(1, z).fdiff()
1211 def _eval_is_real(self):
1212 z = self.args[0]
1213 return polygamma(1, z).is_real
1215 def _eval_is_positive(self):
1216 z = self.args[0]
1217 return polygamma(1, z).is_positive
1219 def _eval_is_negative(self):
1220 z = self.args[0]
1221 return polygamma(1, z).is_negative
1223 def _eval_aseries(self, n, args0, x, logx):
1224 as_polygamma = self.rewrite(polygamma)
1225 args0 = [S.One,] + args0
1226 return as_polygamma._eval_aseries(n, args0, x, logx)
1228 @classmethod
1229 def eval(cls, z):
1230 return polygamma(1, z)
1232 def _eval_expand_func(self, **hints):
1233 z = self.args[0]
1234 return polygamma(1, z).expand(func=True)
1236 def _eval_rewrite_as_zeta(self, z, **kwargs):
1237 return zeta(2, z)
1239 def _eval_rewrite_as_polygamma(self, z, **kwargs):
1240 return polygamma(1, z)
1242 def _eval_rewrite_as_harmonic(self, z, **kwargs):
1243 return -harmonic(z - 1, 2) + pi**2 / 6
1245 def _eval_as_leading_term(self, x, logx=None, cdir=0):
1246 z = self.args[0]
1247 return polygamma(1, z).as_leading_term(x)
1250###############################################################################
1251##################### COMPLETE MULTIVARIATE GAMMA FUNCTION ####################
1252###############################################################################
1255class multigamma(Function):
1256 r"""
1257 The multivariate gamma function is a generalization of the gamma function
1259 .. math::
1260 \Gamma_p(z) = \pi^{p(p-1)/4}\prod_{k=1}^p \Gamma[z + (1 - k)/2].
1262 In a special case, ``multigamma(x, 1) = gamma(x)``.
1264 Examples
1265 ========
1267 >>> from sympy import S, multigamma
1268 >>> from sympy import Symbol
1269 >>> x = Symbol('x')
1270 >>> p = Symbol('p', positive=True, integer=True)
1272 >>> multigamma(x, p)
1273 pi**(p*(p - 1)/4)*Product(gamma(-_k/2 + x + 1/2), (_k, 1, p))
1275 Several special values are known:
1277 >>> multigamma(1, 1)
1278 1
1279 >>> multigamma(4, 1)
1280 6
1281 >>> multigamma(S(3)/2, 1)
1282 sqrt(pi)/2
1284 Writing ``multigamma`` in terms of the ``gamma`` function:
1286 >>> multigamma(x, 1)
1287 gamma(x)
1289 >>> multigamma(x, 2)
1290 sqrt(pi)*gamma(x)*gamma(x - 1/2)
1292 >>> multigamma(x, 3)
1293 pi**(3/2)*gamma(x)*gamma(x - 1)*gamma(x - 1/2)
1295 Parameters
1296 ==========
1298 p : order or dimension of the multivariate gamma function
1300 See Also
1301 ========
1303 gamma, lowergamma, uppergamma, polygamma, loggamma, digamma, trigamma,
1304 beta
1306 References
1307 ==========
1309 .. [1] https://en.wikipedia.org/wiki/Multivariate_gamma_function
1311 """
1312 unbranched = True
1314 def fdiff(self, argindex=2):
1315 from sympy.concrete.summations import Sum
1316 if argindex == 2:
1317 x, p = self.args
1318 k = Dummy("k")
1319 return self.func(x, p)*Sum(polygamma(0, x + (1 - k)/2), (k, 1, p))
1320 else:
1321 raise ArgumentIndexError(self, argindex)
1323 @classmethod
1324 def eval(cls, x, p):
1325 from sympy.concrete.products import Product
1326 if p.is_positive is False or p.is_integer is False:
1327 raise ValueError('Order parameter p must be positive integer.')
1328 k = Dummy("k")
1329 return (pi**(p*(p - 1)/4)*Product(gamma(x + (1 - k)/2),
1330 (k, 1, p))).doit()
1332 def _eval_conjugate(self):
1333 x, p = self.args
1334 return self.func(x.conjugate(), p)
1336 def _eval_is_real(self):
1337 x, p = self.args
1338 y = 2*x
1339 if y.is_integer and (y <= (p - 1)) is True:
1340 return False
1341 if intlike(y) and (y <= (p - 1)):
1342 return False
1343 if y > (p - 1) or y.is_noninteger:
1344 return True