sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.49k stars 4.32k forks source link

[WIP] Fix digamma integral #26563

Open rogerbalsach opened 4 weeks ago

rogerbalsach commented 4 weeks ago

References to other Issues or PRs

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(a0)$$ $$\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 s0 \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) + \gammaE$$ 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

sympy-bot commented 4 weeks ago

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.

:x: 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. #### References to other Issues or PRs 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 * functions.special * Implemented `lerchphi._eval_is_finite` * Implemented `lerchphi._eval_as_leading_term` * Added functionality to `zeta._eval_as_leading_term`