[WIP] Fix digamma integral #26563

rogerbalsach · 2024-05-04T15:48:02Z

References to other Issues or PRs

Brief description of what is fixed or changed

Added regression tests
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:
$$\Phi(0, s, a) = \frac{1}{a^s}$$
Which diverges when $a=0$ and $s \in \mathbb{R}^+$.
When $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}$$
Finally, for $z=1$, the function has the same poles as the zeta function: $s = 1$.
No analytic continuation has been considered, except for $\Phi(1,s,a)=\zeta(s,a)$.
lerchphi._eval_as_leading_term has been implemented. It handles the following limits:
$$\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)$$
Added new functionality to the zeta._eval_as_leading_term method.
It adds the following cases:
$$\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}$$
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$$
We need to correctly compute the following limit:
$$\lim_{t \to 1^-} \left[-t^{z+1}\Phi(t, 1, z+1) - \log(t - 1)\right] = \psi(z+1) + \gamma_E - i \pi$$
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 $\Phi(1, 1, z+1)$ is infinite to avoid simply substituting $t=1$. This requires the implementation of 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 the heuristics algorithm, which also has this bug.
The heuristics algorithm seems to just do the substitution $t = 1$. So a check to make sure the result is finite is needed.
With 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 lerchphi._eval_as_leading_term is needed, which should be able to correctly identify
$$\lim_{\varepsilon \to 0^+} \Phi(1-\varepsilon, 1, z+1) = -\log(\varepsilon) - \gamma_E - \psi(z+1)$$
Finally, a simplification needs to be added somewhere to simplify

(-t*t**z*z*lerchphi(t, 1, z + 1)*gamma(z + 1)/gamma(z + 2)
 - t*t**z*lerchphi(t, 1, z + 1)*gamma(z + 1)/gamma(z + 2)
 - log(t - 1))

to

-t**(z + 1)*lerchphi(t, 1, z + 1) - log(t - 1)

Release Notes

functions.special
- Implemented lerchphi._eval_is_finite
- Implemented lerchphi._eval_as_leading_term
- Added functionality to zeta._eval_as_leading_term

sympy-bot · 2024-05-04T15:48:09Z

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.

functions.special is not a valid release notes header. Release notes headers should be SymPy submodule names, like
```
* core
  * made Add faster

* printing
  * improve LaTeX printing of fractions
```
or other.

If you have added a new submodule, please add it to the list of valid release notes headers at https://github.com/sympy/sympy-bot/blob/master/sympy_bot/submodules.txt.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #26523


#### Brief description of what is fixed or changed

Added regression tests
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:
$$\Phi(0, s, a) = \frac{1}{a^s}$$
Which diverges when $a=0$ and $s \in \mathbb{R}^+$.
When $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}$$
Finally, for $z=1$, the function has the same poles as the zeta function: $s = 1$.
No analytic continuation has been considered, except for $\Phi(1,s,a)=\zeta(s,a)$.
`lerchphi._eval_as_leading_term` has been implemented. It handles the following limits:
$$\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)$$
Added new functionality to the `zeta._eval_as_leading_term` method.
It adds the following cases:
$$\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}$$
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$$
We need to correctly compute the following limit:
$$\lim\_{t \to 1^-} \left\[-t^{z+1}\Phi(t, 1, z+1) - \log(t - 1)\right\] = \psi(z+1) + \gamma_E - i \pi$$
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 $\Phi(1, 1, z+1)$ is infinite to avoid simply substituting $t=1$. This requires the implementation of `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 the `heuristics` algorithm, which also has this bug.
The `heuristics` algorithm seems to just do the substitution $t = 1$. So a check to make sure the result is finite is needed.
With 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 `lerchphi._eval_as_leading_term` is needed, which should be able to correctly identify
$$\lim\_{\varepsilon \to 0^+} \Phi(1-\varepsilon, 1, z+1) = -\log(\varepsilon) - \gamma_E - \psi(z+1)$$
Finally, a simplification needs to be added somewhere to simplify
```python
(-t*t**z*z*lerchphi(t, 1, z + 1)*gamma(z + 1)/gamma(z + 2)
 - t*t**z*lerchphi(t, 1, z + 1)*gamma(z + 1)/gamma(z + 2)
 - log(t - 1))
```
to
```python
-t**(z + 1)*lerchphi(t, 1, z + 1) - log(t - 1)
```

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->

* functions.special
  * Implemented `lerchphi._eval_is_finite`
  * Implemented `lerchphi._eval_as_leading_term`
  * Added functionality to `zeta._eval_as_leading_term`

<!-- END RELEASE NOTES -->

oscarbenjamin · 2024-05-04T18:58:43Z

sympy/functions/special/zeta_functions.py

+    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


What should happen here if any of the is_ checks gives None?

https://docs.sympy.org/latest/guides/booleans.html

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?

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?

It is because a lot of situations where the function is infinite is because one has $0^{-s}$. The previous check 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.

oscarbenjamin · 2024-05-05T13:16:52Z

sympy/functions/special/zeta_functions.py

+        exp_expr = self.expand(func=True)
+        if exp_expr != self:
+            return exp_expr.is_finite


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.

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 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

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?

Yes, the checks need to be from scratch. It is also acceptable that some assumptions queries just return None.

sympy/functions/special/zeta_functions.py

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.

Add regression tests

56290f8

Implement LerchPhi.is_finite

c1db033

oscarbenjamin reviewed May 4, 2024

View reviewed changes

rogerbalsach marked this pull request as ready for review May 4, 2024 22:44

rogerbalsach marked this pull request as draft May 4, 2024 22:45

Fixed bugs with LerchPhi.is_finite

aa567b2

Fixed some edge cases. Check that s is real for 0**s. Return None when the arguments are infinite.

oscarbenjamin reviewed May 5, 2024

View reviewed changes

rogerbalsach commented May 6, 2024

View reviewed changes

sympy/functions/special/zeta_functions.py Show resolved Hide resolved

Remove expand and fix bug.

03693d9

Remove the expand method to avoid a possible infinite recusion. Fixed a bug and added the corresponding test.

oscarbenjamin added assumptions limits labels May 29, 2024

Implemented lerchphi._eval_as_leading_term

85de79f

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.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] Fix digamma integral #26563

[WIP] Fix digamma integral #26563

rogerbalsach commented May 4, 2024 •

edited

sympy-bot commented May 4, 2024 •

edited

oscarbenjamin May 4, 2024

rogerbalsach May 4, 2024

oscarbenjamin May 4, 2024

rogerbalsach May 5, 2024

oscarbenjamin May 5, 2024

rogerbalsach May 6, 2024

oscarbenjamin May 6, 2024

rogerbalsach May 28, 2024

oscarbenjamin May 29, 2024

[WIP] Fix digamma integral #26563

Are you sure you want to change the base?

[WIP] Fix digamma integral #26563

Conversation

rogerbalsach commented May 4, 2024 • edited

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

sympy-bot commented May 4, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rogerbalsach commented May 4, 2024 •

edited

sympy-bot commented May 4, 2024 •

edited