-
-
Notifications
You must be signed in to change notification settings - Fork 4.3k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
[WIP] Fix digamma integral #26563
base: master
Are you sure you want to change the base?
[WIP] Fix digamma integral #26563
Conversation
Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes. ❌ There was an issue with the release notes. Please do not close this pull request; instead edit the description after reading the guide on how to write release notes.
Click here to see the pull request description that was parsed.
|
def _eval_is_finite(self): | ||
z, s, a = self.args | ||
if z == 0 or a == 0: | ||
return z.is_zero and a.is_zero and s.is_nonpositive |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What should happen here if any of the is_
checks gives None
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi @oscarbenjamin , thanks for taking a look and the information about the booleans.
I thought python and was compatible with the fuzzy logic but it is indeed not. I will use the fuzzy_and function, this should handle correctly the None cases.
Also, I am not sure what is the convention when the function evaluates to NaN, should the is_finite be None?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
what is the convention when the function evaluates to NaN, should the is_finite be None?
Probably. What case evaluates to nan that you are thinking of?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It is because a lot of situations where the function is infinite is because one has s.is_nonpositive
gives False when s is not real. I have added a check for s to be real, when s is complex the expression is undefined so now the function should return None.
Fixed some edge cases. Check that s is real for 0**s. Return None when the arguments are infinite.
exp_expr = self.expand(func=True) | ||
if exp_expr != self: | ||
return exp_expr.is_finite |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What is expand(func=True)
being used for here?
We should avoid creating new expressions during an assumptions query. Just creating the new expression runs the risk of triggering an identical assumptions query which leads to infinite recursion.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, I did not think about that. My idea was that, since LerchPhi includes a lot of common functions as particular values (zeta function, polylog, Dirichlet eta, ...), I can try to write the function in terms of these other functions and then let them determine whether the function is finite or not. Is there a safe way to do it without the risks you are mentioning?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There is not really a safe way to do it. We just have to accept that the core assumptions cannot resolve all queries:
https://docs.sympy.org/latest/guides/assumptions.html
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, sorry for the delay, but I haven't had time these last few weeks. I understand then that this part of the code should be changed and all the checks need to be implemented from scratch. Is that correct?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, the checks need to be from scratch. It is also acceptable that some assumptions queries just return None.
Remove the expand method to avoid a possible infinite recusion. Fixed a bug and added the corresponding test.
Implemented functions lerchphi._eval_as_leading_term for a few selected cases. Added new cases to the function zeta._eval_as_leading_term needed for the lerchphi implementation. Added tests for both implementations.
References to other Issues or PRs
Fixes #26523
Brief description of what is fixed or changed
Added regression tests
$$\Phi(0, s, a) = \frac{1}{a^s}$$ $a=0$ and $s \in \mathbb{R}^+$ .$z \neq 0$ , the function diverges when $a \in - \mathbb{N}$ and $s \in \mathbb{R}^+$ .
$$\Phi(z, s, -k) = \sum_{n=0}^\infty \frac{z^n}{(n-k)^s}$$ $z=1$ , the function has the same poles as the zeta function: $s = 1$ .$\Phi(1,s,a)=\zeta(s,a)$ .
$$\lim_{z \to 1} \Phi(z, 1, a_0) = -\log(\log(z)) - \gamma_E - \psi(a_0)$$
$$\lim_{z \to 1} \Phi(z, s, a) = \zeta(s, a)$$
$$\lim_{(s,a) \to (s_0, a_0)} \zeta(s, a) = \zeta(s_0, a_0), \qquad \text{if} \qquad s_0 \neq 1$$
$$\lim_{(s,a) \to (1, a_0)} \zeta(s, a) = \frac{1}{s-1}$$
Implemented the functions _eval_is_finite and _eval_as_leading_term for the lerchphi class.
lerchphi._eval_is_finite
has been implemented with finite complex numbers in mind. A possible extension to handle infinite arguments is left for the future.The divergences considered are the following:
Which diverges when
When
Finally, for
No analytic continuation has been considered, except for
lerchphi._eval_as_leading_term
has been implemented. It handles the following limits:Added new functionality to the
zeta._eval_as_leading_term
method.It adds the following cases:
Added tests for all the new changes.
Other comments
Road map:
Following the discussion at https://github.com/sympy/sympy/issues/26523, in order to solve the bug with the integral
$$\int_0^1 \frac{1-t^z}{1-t} \mathrm{d}t = \psi(z+1) + \gamma_E$$
$$\lim_{t \to 1^-} \left[-t^{z+1}\Phi(t, 1, z+1) - \log(t - 1)\right] = \psi(z+1) + \gamma_E - i \pi$$ $\Phi(1, 1, z+1)$ is infinite to avoid simply substituting $t=1$ . This requires the implementation of $t = 1$ . So a check to make sure the result is finite is needed.
We need to correctly compute the following limit:
Sympy currently gives the wrong result for this limit, which causes the bug. To avoid returning the wrong value, we need the following:
We first need sympy to identify that
lerchphi._eval_is_finite
.This is enough to solve the bug in the
gruntz
algorithm. However, Sympy will also try to solve the limit using theheuristics
algorithm, which also has this bug.The
heuristics
algorithm seems to just do the substitutionWith these, the integral should no longer return the wrong result but rather either raise an error or return the unevaluated result.
In order to compute the correct limit, an implementation of
$$\lim_{\varepsilon \to 0^+} \Phi(1-\varepsilon, 1, z+1) = -\log(\varepsilon) - \gamma_E - \psi(z+1)$$
lerchphi._eval_as_leading_term
is needed, which should be able to correctly identifyFinally, a simplification needs to be added somewhere to simplify
to
Release Notes
lerchphi._eval_is_finite
lerchphi._eval_as_leading_term
zeta._eval_as_leading_term