Implement digamma and trigamma functions as Function subclasses #17615
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
✅ Hi, I am the SymPy bot (v149). I'm here to help you write a release notes entry. Please read the guide on how to write release notes. Your release notes are in good order. Here is what the release notes will look like: This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.5. Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. Click here to see the pull request description that was parsed.
Update The release notes on the wiki have been updated. |
czgdp1807
reviewed
Sep 15, 2019
czgdp1807
reviewed
Sep 15, 2019
czgdp1807
reviewed
Sep 15, 2019
czgdp1807
reviewed
Sep 15, 2019
czgdp1807
reviewed
Sep 15, 2019
|
@asmeurer Please review this PR. I have provided some of the comments from my side. |
|
This may fix the problems with gruntz test. $ git diff
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index d98b3cd006..3d9d2ab725 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -1048,43 +1048,9 @@ def _eval_is_negative(self):
def _eval_aseries(self, n, args0, x, logx):
from sympy import Order
- if args0[1] != oo:
- return super(polygamma, self)._eval_aseries(n, args0, x, logx)
- z = self.args[0]
- N = sympify(0)
-
- if N == 0:
- # digamma function series
- # Abramowitz & Stegun, p. 259, 6.3.18
- r = log(z) - 1/(2*z)
- o = None
- if n < 2:
- o = Order(1/z, x)
- else:
- m = ceiling((n + 1)//2)
- l = [bernoulli(2*k) / (2*k*z**(2*k)) for k in range(1, m)]
- r -= Add(*l)
- o = Order(1/z**(2*m), x)
- return r._eval_nseries(x, n, logx) + o
- else:
- # proper polygamma function
- # Abramowitz & Stegun, p. 260, 6.4.10
- # We return terms to order higher than O(x**n) on purpose
- # -- otherwise we would not be able to return any terms for
- # quite a long time!
- fac = gamma(N)
- e0 = fac + N*fac/(2*z)
- m = ceiling((n + 1)//2)
- for k in range(1, m):
- fac = fac*(2*k + N - 1)*(2*k + N - 2) / ((2*k)*(2*k - 1))
- e0 += bernoulli(2*k)*fac/z**(2*k)
- o = Order(1/z**(2*m), x)
- if n == 0:
- o = Order(1/z, x)
- elif n == 1:
- o = Order(1/z**2, x)
- r = e0._eval_nseries(z, n, logx) + o
- return (-1 * (-1/z)**N * r)._eval_nseries(x, n, logx)
+ as_polygamma = self.rewrite(polygamma)
+ args0 = [S.Zero,] + args0
+ return as_polygamma._eval_aseries(n, args0, x, logx)
@classmethod
def eval(cls, z):
@@ -1189,6 +1155,9 @@ def _eval_expand_func(self, **hints):
def _eval_rewrite_as_harmonic(self, z, **kwargs):
return harmonic(z - 1) - S.EulerGamma
+ def _eval_rewrite_as_polygamma(self, z, **kwargs):
+ return polygamma(0, z)
+
def _eval_as_leading_term(self, x):
from sympy import Order
n = sympify(0)I think that it would still be a better idea to make |
|
@sylee957 I was not able to call polygamma's private methods from inside digamma or trigamma. That's why i had to duplicate their implementation. |
Codecov Report
@@ Coverage Diff @@
## master #17615 +/- ##
============================================
- Coverage 74.803% 74.68% -0.124%
============================================
Files 633 634 +1
Lines 164986 165660 +674
Branches 38791 38961 +170
============================================
+ Hits 123416 123716 +300
- Misses 36169 36453 +284
- Partials 5401 5491 +90 |
Code
$ git diff
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index 37df26ee78..b1c03756f9 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -1020,153 +1020,9 @@ class digamma(Function):
.. [2] http://mathworld.wolfram.com/DigammaFunction.html
.. [3] http://functions.wolfram.com/GammaBetaErf/PolyGamma2/
"""
- def _eval_evalf(self, prec):
- n = sympify(0)
- # the mpmath polygamma implementation valid only for nonnegative integers
- if n.is_number and n.is_real:
- if (n.is_integer or n == int(n)) and n.is_nonnegative:
- return super(digamma, self)._eval_evalf(prec)
-
- def fdiff(self, argindex=1):
- if argindex == 1:
- x = self.args[0]
- return trigamma(x)
- else:
- raise ArgumentIndexError(self, argindex)
-
- def _eval_is_real(self):
- if self.args[0].is_positive:
- return True
-
- def _eval_is_positive(self):
- if self.args[0].is_positive:
- return False
-
- def _eval_is_negative(self):
- if self.args[0].is_positive:
- return True
-
- def _eval_aseries(self, n, args0, x, logx):
- from sympy import Order
- as_polygamma = self.rewrite(polygamma)
- args0 = [sympify(0),] + args0
- return as_polygamma._eval_aseries(n,args0,x,logx)
-
@classmethod
def eval(cls, z):
- n, z = map(sympify, (0, z))
- from sympy import unpolarify
-
-
- if n.is_integer:
- if n.is_nonnegative:
- nz = unpolarify(z)
- if z != nz:
- return polygamma(n, nz)
-
- if n == -1:
- return loggamma(z)
- else:
- if z.is_Number:
- if z is S.NaN:
- return S.NaN
- elif z is S.Infinity:
- if n.is_Number:
- if n is S.Zero:
- return S.Infinity
- else:
- return S.Zero
- elif z.is_Integer:
- if z.is_nonpositive:
- return S.ComplexInfinity
- else:
- if n is S.Zero:
- return -S.EulerGamma + harmonic(z - 1, 1)
- elif n.is_odd:
- return (-1)**(n + 1)*factorial(n)*zeta(n + 1, z)
-
- if n == 0:
- if z is S.NaN:
- return S.NaN
- elif z.is_Rational:
-
- p, q = z.as_numer_denom()
-
- # only expand for small denominators to avoid creating long expressions
- if q <= 5:
- return expand_func(polygamma(n, z, evaluate=False))
-
- elif z in (S.Infinity, S.NegativeInfinity):
- return S.Infinity
- else:
- t = z.extract_multiplicatively(S.ImaginaryUnit)
- if t in (S.Infinity, S.NegativeInfinity):
- return S.Infinity
-
- # TODO n == 1 also can do some rational z
-
- def _eval_expand_func(self, **hints):
- n = sympify(0)
- z = self.args[0]
-
- if n.is_Integer and n.is_nonnegative:
- if z.is_Add:
- coeff = z.args[0]
- if coeff.is_Integer:
- e = -(n + 1)
- if coeff > 0:
- tail = Add(*[Pow(
- z - i, e) for i in range(1, int(coeff) + 1)])
- else:
- tail = -Add(*[Pow(
- z + i, e) for i in range(0, int(-coeff))])
- return polygamma(n, z - coeff) + (-1)**n*factorial(n)*tail
-
- elif z.is_Mul:
- coeff, z = z.as_two_terms()
- if coeff.is_Integer and coeff.is_positive:
- tail = [ polygamma(n, z + Rational(
- i, coeff)) for i in range(0, int(coeff)) ]
- if n == 0:
- return Add(*tail)/coeff + log(coeff)
- else:
- return Add(*tail)/coeff**(n + 1)
- z *= coeff
-
- if n == 0 and z.is_Rational:
- p, q = z.as_numer_denom()
-
- # Reference:
- # Values of the polygamma functions at rational arguments, J. Choi, 2007
- part_1 = -S.EulerGamma - pi * cot(p * pi / q) / 2 - log(q) + Add(
- *[cos(2 * k * pi * p / q) * log(2 * sin(k * pi / q)) for k in range(1, q)])
-
- if z > 0:
- n = floor(z)
- z0 = z - n
- return part_1 + Add(*[1 / (z0 + k) for k in range(n)])
- elif z < 0:
- n = floor(1 - z)
- z0 = z + n
- return part_1 - Add(*[1 / (z0 - 1 - k) for k in range(n)])
-
- return polygamma(n, z)
-
- def _eval_rewrite_as_harmonic(self, z, **kwargs):
- return harmonic(z - 1) - S.EulerGamma
-
- def _eval_rewrite_as_polygamma(self, z, **kwargs):
- return polygamma(0,z)
-
- def _eval_as_leading_term(self, x):
- from sympy import Order
- n = sympify(0)
- z = self.args[0].as_leading_term(x)
- o = Order(z, x)
- if n == 0 and o.contains(1/x):
- return o.getn() * log(x)
- else:
- return self.func(n, z)
+ return polygamma(0, z)
class trigamma(Function):
@@ -1196,152 +1052,9 @@ class trigamma(Function):
.. [2] http://mathworld.wolfram.com/TrigammaFunction.html
.. [3] http://functions.wolfram.com/GammaBetaErf/PolyGamma2/
"""
- def _eval_evalf(self, prec):
- n = sympify(1)
- # the mpmath polygamma implementation valid only for nonnegative integers
- if n.is_number and n.is_real:
- if (n.is_integer or n == int(n)) and n.is_nonnegative:
- return super(trigamma, self)._eval_evalf(prec)
-
- def fdiff(self, argindex=1):
- if argindex == 1:
- n = sympify(1)
- z = self.args[0]
- return polygamma(n + 1, z)
- else:
- raise ArgumentIndexError(self, argindex)
-
- def _eval_is_real(self):
- if self.args[0].is_positive:
- return True
-
- def _eval_is_positive(self):
- if self.args[0].is_positive:
- return True
-
- def _eval_is_negative(self):
- if self.args[0].is_positive:
- return False
-
- def _eval_aseries(self, n, args0, x, logx):
- from sympy import Order
- as_polygamma = self.rewrite(polygamma)
- args0 = [S.One,] + args0
- return as_polygamma._eval_aseries(n,args0,x,logx)
-
@classmethod
def eval(cls, z):
- n, z = map(sympify, (1, z))
- from sympy import unpolarify
-
- if n.is_integer:
- if n.is_nonnegative:
- nz = unpolarify(z)
- if z != nz:
- return polygamma(n, nz)
-
- if n == -1:
- return loggamma(z)
- else:
- if z.is_Number:
- if z is S.NaN:
- return S.NaN
- elif z is S.Infinity:
- if n.is_Number:
- if n is S.Zero:
- return S.Infinity
- else:
- return S.Zero
- elif z.is_Integer:
- if z.is_nonpositive:
- return S.ComplexInfinity
- else:
- if n is S.Zero:
- return -S.EulerGamma + harmonic(z - 1, 1)
- elif n.is_odd:
- return (-1)**(n + 1)*factorial(n)*zeta(n + 1, z)
-
- if n == 0:
- if z is S.NaN:
- return S.NaN
- elif z.is_Rational:
-
- p, q = z.as_numer_denom()
-
- # only expand for small denominators to avoid creating long expressions
- if q <= 5:
- return expand_func(polygamma(n, z, evaluate=False))
-
- elif z in (S.Infinity, S.NegativeInfinity):
- return S.Infinity
- else:
- t = z.extract_multiplicatively(S.ImaginaryUnit)
- if t in (S.Infinity, S.NegativeInfinity):
- return S.Infinity
-
- # TODO n == 1 also can do some rational z
-
- def _eval_expand_func(self, **hints):
- n = sympify(1)
- z = self.args[0]
-
- if n.is_Integer and n.is_nonnegative:
- if z.is_Add:
- coeff = z.args[0]
- if coeff.is_Integer:
- e = -(n + 1)
- if coeff > 0:
- tail = Add(*[Pow(
- z - i, e) for i in range(1, int(coeff) + 1)])
- else:
- tail = -Add(*[Pow(
- z + i, e) for i in range(0, int(-coeff))])
- return polygamma(n, z - coeff) + (-1)**n*factorial(n)*tail
-
- elif z.is_Mul:
- coeff, z = z.as_two_terms()
- if coeff.is_Integer and coeff.is_positive:
- tail = [ polygamma(n, z + Rational(
- i, coeff)) for i in range(0, int(coeff)) ]
- if n == 0:
- return Add(*tail)/coeff + log(coeff)
- else:
- return Add(*tail)/coeff**(n + 1)
- z *= coeff
-
- if n == 0 and z.is_Rational:
- p, q = z.as_numer_denom()
-
- # Reference:
- # Values of the polygamma functions at rational arguments, J. Choi, 2007
- part_1 = -S.EulerGamma - pi * cot(p * pi / q) / 2 - log(q) + Add(
- *[cos(2 * k * pi * p / q) * log(2 * sin(k * pi / q)) for k in range(1, q)])
-
- if z > 0:
- n = floor(z)
- z0 = z - n
- return part_1 + Add(*[1 / (z0 + k) for k in range(n)])
- elif z < 0:
- n = floor(1 - z)
- z0 = z + n
- return part_1 - Add(*[1 / (z0 - 1 - k) for k in range(n)])
-
- return polygamma(n, z)
-
- def _eval_rewrite_as_zeta(self, z, **kwargs):
- return (-1)**(2)*factorial(1)*zeta(2, z)
-
- def _eval_rewrite_as_polygamma(self, z, **kwargs):
- return polygamma(0,z)
-
- def _eval_rewrite_as_harmonic(self, n, z, **kwargs):
- return S.NegativeOne**(2) * factorial(1) * (zeta(2) - harmonic(0, 2))
-
- def _eval_as_leading_term(self, x):
- from sympy import Order
- n = sympify(1)
- z = self.args[0].as_leading_term(x)
- return self.func(n, z)
+ return polygamma(1, z)Actually, a lot of things to be easier if these trivial gamma functions can automatically canonicalize as |
|
@sylee957 is this what you were expecting? |
sylee957
reviewed
Sep 20, 2019
|
|
||
| def _eval_is_real(self): | ||
| z = self.args[0] | ||
| return polygamma(0,z).is_real() |
There was a problem hiding this comment.
This should be polygamma(1, z)
sylee957
reviewed
Sep 20, 2019
There was a problem hiding this comment.
I don't think that it is necessary to implement every computation methods that polygamma have, if we're only going to make digamma and trigamma as notational purposes for a while.
Of course, adding these methods can be also be fine if you want to keep the consistency of assumptions or such, but one thing you have to consider is that you have to add the test cases for all these methods to make sure that they work properly.
You can add tests like digamma(z, evaluate=False).is_positive or such.
I also have some doubts that digamma(z, evaluate=False).is_real() can work for now.
|
@sylee957 Yes, But it doesn't seem to be working Any suggestions? |
|
@sylee957 Please address #17615 (comment) |
|
You may have to delete the |
|
@sylee957 Is this okay? |
sylee957
reviewed
Oct 1, 2019
sylee957
reviewed
Oct 1, 2019
Co-Authored-By: S.Y. Lee <sylee957@gmail.com>
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
References to other Issues or PRs
Fixes #17596
Brief description of what is fixed or changed
The digamma and trigamma functions were Python functions. Now they are Function subclasses
Other comments
Added tests for digamma and trigamma
Release Notes