B(1) = +½ set 2: generalised Bernoulli function

[Parcly-Taxel]

References to other Issues or PRs

Part 2 of the changes formerly included in #23926, building upon #23973. The changes here are also the actual changes that will fix #23866.

Brief description of what is fixed or changed

This PR generalises bernoulli() to the generalised Bernoulli function described in Peter Luschny's "An introduction to the Bernoulli function":
$$B(s,a)=\begin{cases}1&s=0\\-s\zeta(1-s,a)&s\ne0\end{cases}$$
In the process it also changes bernoulli(1) to +1/2.

Other comments

Since the Genocchi numbers are currently defined in terms of the Bernoulli numbers, this also changes genocchi(1) to -1; a later PR will generalise this integer sequence to the two-complex-argument generalised Genocchi function in the same vein as the generalised Bernoulli function, though.

Luschny also defines the central Bernoulli numbers and polynomials, but I have not made a new SymPy function out of that because its generalisation to arbitrary complex arguments can be expressed as $2^s B(s,\frac{a+1}2)$.

Release Notes

• functions
  • BREAKING: bernoulli(1) is now 1/2 instead of -1/2 and genocchi(1) is now -1 instead of 1.
  • bernoulli() now accepts complex arguments, implementing Luschny (2020).

[sympy-bot]

✅

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:

• functions

  • BREAKING: bernoulli(1) is now 1/2 instead of -1/2 and genocchi(1) is now -1 instead of 1. (#23984 by @Parcly-Taxel)

  • bernoulli() now accepts complex arguments, implementing Luschny (2020). (#23984 by @Parcly-Taxel)

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.

#### References to other Issues or PRs

Part 2 of the changes formerly included in #23926, building upon #23973. The changes here are also the actual changes that will fix #23866.

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

This PR generalises `bernoulli()` to the generalised Bernoulli function described in Peter Luschny's "[An introduction to the Bernoulli function](https://arxiv.org/abs/2009.06743)":
$$B(s,a)=\begin{cases}1&s=0\\\\-s\zeta(1-s,a)&s\ne0\end{cases}$$
In the process it also changes `bernoulli(1)` to `+1/2`.

#### Other comments

Since the Genocchi numbers are currently defined in terms of the Bernoulli numbers, this also changes `genocchi(1)` to `-1`; a later PR will generalise this integer sequence to the two-complex-argument generalised Genocchi function in the same vein as the generalised Bernoulli function, though.

Luschny also defines the _central Bernoulli numbers_ and _polynomials_, but I have not made a new SymPy function out of that because its generalisation to arbitrary complex arguments can be expressed as $2^s B(s,\frac{a+1}2)$.

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* functions
    * BREAKING: `bernoulli(1)` is now `1/2` instead of `-1/2` and `genocchi(1)` is now `-1` instead of `1`.
    * `bernoulli()` now accepts complex arguments, implementing [Luschny (2020)](https://arxiv.org/abs/2009.06743).
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

       before           after         ratio
     [41d90958]       [a3826118]
-         987±3μs          626±2μs     0.63  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-     2.84±0.01ms         1.16±0ms     0.41  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     5.62±0.01ms         1.70±0ms     0.30  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[oscarbenjamin]

I think this looks good. Does anyone else have any opinions about this?

The discussion for the rationale for this change is in the linked issue and reference cited above.

Not a new problem but this should also get fixed:

In [1]: lambdify(x, bernoulli(x))(1.2)
Out[1]: array([ 1. , -0.5])

That's scipy's bernoulli function which is inconsistent with SymPy's.

[Parcly-Taxel]

Not a new problem but this should also get fixed: … That's scipy's bernoulli function which is inconsistent with SymPy's.

I think I've got the solution for that. Put the following code in sympy/printing/numpy.py:

    def _print_bernoulli(self, expr):
        return self._print(expr._eval_rewrite_as_zeta(*expr.args))

Does this make sense and should I commit it in this PR?

[oscargus]

I think that makes sense and could go into this. Possibly with a small comment of the type # scipy's bernoulli in inconsistent with SymPy's so rewrite and a test. Thanks!

[oscarbenjamin]

Okay let's get this in.

Thanks @Parcly-Taxel!

[asmeurer]

The lambdified bernoulli was already wrong anyway. scipy.special.bernoulli(n) returns the first n Bernoulli numbers. The equivalent for SymPy would be scipy.special.bernoulli(n)[-1].

[fchapoton]

This is a very wrong move and should be undone. You have been tricked to do that by people with no authority whatsoever. This is not good to try changing the stable convention used by the large majority of mathematical papers.

[Parcly-Taxel]

This is a very wrong move and should be undone. You have been tricked to do that by people with no authority whatsoever. This is not good to try changing the stable convention used by the large majority of mathematical papers.

You have no authority here to decide what $B_1$ is. Meurer, Benjamin and I have already decided that it will change to $+\frac12$ in 1.12 and that such change is entirely OK; I layered my PRs in such a way to lock this change in. It's too late to turn back.

[edgarcosta]

@Parcly-Taxel, I think you have just agreed that you tricked them:

I layered my PRs in such a way to lock this change in. It's too late to turn back.

[oscarbenjamin]

I layered my PRs in such a way to lock this change in. It's too late to turn back.

It's not too late to revert anything. Nothing has been released yet.

[Parcly-Taxel]

I think you have just agreed that you tricked them:

There was no trickery involved here. There is also no authority defining what the Bernoulli numbers are, just as there is no authority defining $\pi$ or the real numbers. I merely wrote code illuminating mathematical truth, which every CAS should strive to do.

[asmeurer]

Can we please keep the discussions here civil and within our code of conduct. Keep the discussion to the technical merits and specific mathematical details. Claiming that people are "playing tricks", especially without any arguments to back it up, is not a helpful avenue of discussion and is unnecessarily antagonistic.

Also, it would be best to have this discussion at #23866. There's a lot of PRs relating to this issue in SymPy, but that's the main issue for bernoulli(1) and related changes. And as I requested elsewhere, if there are relevant discussions happening on other forums (such as Sage forums), please cross-link them there as many of us are not subscribed to them.

As Oscar stated, this is not yet "locked in", because it hasn't been in a release. That was done on purpose because we knew this could be controversial.