B(1) = +½ set 1: analytically correct Hurwitz zeta #23973
Conversation
In Mathematica there are two two-argument extensions of the Riemann zeta
function, Zeta[s,a] and HurwitzZeta[s,a]. The latter function provides
the correct analytic continuation of the Dirichlet sum to the left
half-plane in a, and thus satisfies a key identity in Apostol's
Introduction to Analytic Number Theory (Theorem 12.13 in the 1995
edition): when n is a non-negative integer, for *any* complex a
-B (a)
n + 1
ζ(-n, a) = ───────────
n + 1
where B denotes the Bernoulli polynomial, e.g. ζ(-1, -4) = -121/12.
There are other nice properties of HurwitzZeta[] which the other
extension Zeta[] lacks. Indeed, Zeta[] mangles the Dirichlet sum on the
left half-plane in a to the point where the graph looks more like
EllipticF[] than Log[].
SymPy's dependency mpmath agrees with HurwitzZeta[], but SymPy itself
has relations implementing Zeta[], so that zeta(-1, -4) = 119/12. This
commit rectifies the discrepancy.
With this change in behaviour also comes a refactoring of the evaluation
tree: ζ(s, a) can be analytically continued to all s ≠ 1 for any a not a
nonpositive integer, while Apostol's identity gives meaningful values to
the holes when s *is* a nonpositive integer. Thus ζ(s, a) is undefined
iff a is a nonpositive integer and s is not. The new evaluation tree
also suppresses explicit a = 1; to effect this change bernoulli()'s
evaluation tree had to be refactored in the same way, but its behaviour
is for the time being unchanged.
|
✅ Hi, I am the SymPy bot (v167). 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.12. Click here to see the pull request description that was parsed.
Update The release notes on the wiki have been updated. |
|
Benchmark results from GitHub Actions Lower numbers are good, higher numbers are bad. A ratio less than 1 Significantly changed benchmark results (PR vs master) Significantly changed benchmark results (master vs previous release) before after ratio
[26f7bdbe] [8bc57b45]
<sympy-1.11^0>
- 963±4μs 614±0.7μs 0.64 solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
- 2.82±0.01ms 1.15±0ms 0.41 solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
- 5.63±0.02ms 1.69±0ms 0.30 solve.TimeSparseSystem.time_linear_eq_to_matrix(30)
Full benchmark results can be found as artifacts in GitHub Actions |
Does anyone else have any opinions on this? I think this change looks good to merge. |
There is, however, a nagging question. Later on I am going to extend This will be used for if and only if However, if you think this PR is good to merge as-is, then merge it. |
These kinds of automatic rewrites are often problematic. Actually I now see that this PR extends this to non-integer In [1]: zeta(S.Half, 2)
Out[1]: -harmonic(1, 1/2) + ζ(1/2)I think it would be better at least not to extend this automatic rewrite. Using this rewrite when both the zeta and the harmonic would fully evaluate so that we don't have any zeta or harmonic functions in the final result probably makes sense. Otherwise I don't think it's a good idea when the arguments are symbolic because of the recursive situations that you describe and also because automatic simplification should only apply when it is unambiguously a "good" simplification. |
In order to prevent potential infinite recursion when harmonic() is extended to complex arguments, the scope of zeta(s,a)'s automatic evaluation when a is an integer has to be restricted. Here the rewrite is performed if and only if s is an Integer and a is a positive Integer, ensuring that the harmonic part collapses into a Rational.
I've sussed out that issue. I have decided that the automatic harmonic rewrite will be applied only if |
|
Okay, I think this is good. Let's get it in. Thanks @Parcly-Taxel! |
In Peter Luschny's "An introduction to the Bernoulli function" (https://arxiv.org/abs/2009.06743, will just be called Luschny in future commits) a generalised (two-argument) Bernoulli function is defined. One of the given representations is ⎧ 1 for s = 0 B(s, a) = ⎨ ⎩-s⋅ζ(1 - s, a) otherwise When s is a nonnegative integer this function reproduces the Bernoulli polynomials; when in addition a = 1 (the default choice for the Hurwitz zeta function's second argument, yielding Riemann's zeta) the Bernoulli numbers are obtained. Hence B(1) = +1/2, which is how Jakob Bernoulli originally defined his numbers in Ars Conjectandi (1713); this choice gives the simpler form/wider range of validity than B(1) = -1/2 for very many equations like Faulhaber's formula and the Euler–Maclaurin summation formula. We use the zeta representation of the Bernoulli function, despite Luschny striving to avoid such dependence by defining it as an integral, since mpmath provides a correct implementation of the Hurwitz zeta function; SymPy's zeta() was also aligned to Hurwitz zeta in sympy#23973.
|
Yes, the best pattern to follow is to only automatically evaluate on explicit Rational inputs, and leave more general symbolic evaluations to other simplification routines like
I'm curious if @jksuom has any thoughts on this and the underlying issue #23866 |
|
I would like to understand why the sign was changed in the first place but otherwise I see no problems with this. |
|
I think the sign was always -1/2 in SymPy. If you mean why it differed historically see the Luschny page linked in the issue. |
References to other Issues or PRs
This is part 1 of the changes formerly included in #23926.
Brief description of what is fixed or changed
mpmath's implementation of the Hurwitz zeta function agrees with Mathematica's
HurwitzZeta[], but the symbolic code in SymPy currently parallelsZeta[]instead. This PR rectifies the discrepancy, especially sinceHurwitzZeta[]is the correct analytic continuation of the Riemann zeta function and thus has the nicer properties.Other comments
bernoulli()is also refactored in this PR, but its user-facing behaviour is not changed;bernoulli(1)still returns-1/2.Release Notes
zeta(s, a)now implements the analytically continued Hurwitz zeta function (Mathematica'sHurwitzZeta[]), not Mathematica'sZeta[].zeta()now simplifies for more sets of parameters.