Repository: sympy/sympy

Issue:
The value of bernoulli(1)
There are two conventions for the Bernoulli numbers, specifically $B_1$ can be defined to be 1/2 or -1/2 (see https://en.wikipedia.org/wiki/Bernoulli_number), depending on whether it is defined using the exponential generating function $\frac{x}{e^x-1}$ or $\frac{xe^x}{e^x - 1}$. SymPy uses the -1/2 definition.

Donald Knuth has recently written that he's been convinced that the -1/2 definition is wrong and historically inaccurate, and he updated *Concrete Mathematics* accordingly (https://www-cs-faculty.stanford.edu/~knuth/news22.html):

> The sequence of Bernoulli numbers plays a prominent role in mathematics, especially in connection with asymptotic expansions related to Euler's summation formula. When I learned about this fascinating sequence, during my undergraduate days, I was taught that $B_1$ is equal to minus one-half. So I duly taught the same to my students, and went on to write books that explained what I thought Bernoulli and Euler had done.
>
> But last year I took a close look at Peter Luschny's [Bernoulli manifesto](http://luschny.de/math/zeta/The-Bernoulli-Manifesto.html), where he gives more than a dozen good reasons why the value of $B_1$ should really be plus one-half. He explains that some mathematicians of the early 20th century had unilaterally changed the conventions, because some of their formulas came out a bit nicer when the negative value was used. It was their well-intentioned but ultimately poor choice that had led to what I'd been taught in the 1950s.
>
> Luschny's webpage cites, for example, recent treatments of the subject by leading mathematicians such as Terence Tao. And his most compelling argument, from my personal perspective, is the way he unveils the early publications: I learned from him that my own presentation of the story, in [The Art of Computer Programming](https://www-cs-faculty.stanford.edu/~knuth/taocp.html) and much more extensively in [Concrete Mathematics](https://www-cs-faculty.stanford.edu/~knuth/gkp.html), was a violation of history! I had put words and thoughts into Bernoulli and Euler's minds that were not theirs at all. This hurt, because I've always tried to present the evolution of ideas faithfully; in this case I'd fooled myself, by trying to conform what they wrote to what I'd learned.
>
> By now, hundreds of books that use the “minus-one-half” convention have unfortunately been written. Even worse, all the major software systems for symbolic mathematics have that 20th-century aberration deeply embedded. Yet Luschny convinced me that we have all been wrong, and that it's high time to change back to the correct definition before the situation gets even worse.
>
> Therefore I changed the definition of $B_1$ in all printings of The Art of Computer Programming during the latter half of 2021. And the new (34th) printing of Concrete Mathematics, released in January 2022, contains the much more extensive changes that are needed to tell a more comprehensive story.

Here is Luschny's page he cites http://luschny.de/math/zeta/The-Bernoulli-Manifesto.html.

So maybe we should reconsider SymPy's definition of `bernoulli`. At the very least, I think there should be a way to specify which convention is used, $B^+_n$ or $B^-_n$.

As an aside, [the docs for `bernoulli`](https://docs.sympy.org/dev/modules/functions/combinatorial.html#sympy.functions.combinatorial.numbers.bernoulli) state that the defining recurrence relation is $$0 = \sum_{k=0}^n \binom{n+1}{k} B_k$$ but for the -1/2 definition it should actually be $$\delta_{m,0} = \sum_{k=0}^n \binom{n+1}{k} B_k$$ (for the +1/2 definition it's $$n + 1 = \sum_{k=0}^n \binom{n+1}{k} B_k,$$ see https://en.wikipedia.org/wiki/Bernoulli_number#Recursive_definition).

Discussion:

Patch:
diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index adf5c56d5c..2a9810a4d1 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -21,7 +21,6 @@
 from mpmath.libmp.backend import MPZ
 from mpmath.libmp.libmpc import _infs_nan
 from mpmath.libmp.libmpf import dps_to_prec, prec_to_dps
-from mpmath.libmp.gammazeta import mpf_bernoulli
 
 from .sympify import sympify
 from .singleton import S
@@ -43,7 +42,6 @@
     from sympy.functions.elementary.complexes import Abs, re, im
     from sympy.functions.elementary.integers import ceiling, floor
     from sympy.functions.elementary.trigonometric import atan
-    from sympy.functions.combinatorial.numbers import bernoulli
     from .numbers import Float, Rational, Integer, AlgebraicNumber, Number
 
 LG10 = math.log(10, 2)
@@ -1046,17 +1044,6 @@ def evalf_piecewise(expr: 'Expr', prec: int, options: OPT_DICT) -> TMP_RES:
     raise NotImplementedError
 
 
-def evalf_bernoulli(expr: 'bernoulli', prec: int, options: OPT_DICT) -> TMP_RES:
-    arg = expr.args[0]
-    if not arg.is_Integer:
-        raise ValueError("Bernoulli number index must be an integer")
-    n = int(arg)
-    b = mpf_bernoulli(n, prec, rnd)
-    if b == fzero:
-        return None, None, None, None
-    return b, None, prec, None
-
-
 def evalf_alg_num(a: 'AlgebraicNumber', prec: int, options: OPT_DICT) -> TMP_RES:
     return evalf(a.to_root(), prec, options)
 
@@ -1400,7 +1387,6 @@ def evalf_symbol(x: 'Expr', prec: int, options: OPT_DICT) -> TMP_RES:
 
 def _create_evalf_table():
     global evalf_table
-    from sympy.functions.combinatorial.numbers import bernoulli
     from sympy.concrete.products import Product
     from sympy.concrete.summations import Sum
     from .add import Add
@@ -1454,7 +1440,6 @@ def _create_evalf_table():
         Product: evalf_prod,
         Piecewise: evalf_piecewise,
 
-        bernoulli: evalf_bernoulli,
         AlgebraicNumber: evalf_alg_num,
     }
 
diff --git a/sympy/functions/combinatorial/numbers.py b/sympy/functions/combinatorial/numbers.py
index e81f5d9413..2b8fe72325 100644
--- a/sympy/functions/combinatorial/numbers.py
+++ b/sympy/functions/combinatorial/numbers.py
@@ -419,23 +419,38 @@ def _eval_rewrite_as_TribonacciConstant(self, n, **kwargs):
 
 class bernoulli(Function):
     r"""
-    Bernoulli numbers / Bernoulli polynomials
+    Bernoulli numbers / Bernoulli polynomials / Bernoulli function
 
     The Bernoulli numbers are a sequence of rational numbers
     defined by `B_0 = 1` and the recursive relation (`n > 0`):
 
-    .. math :: 0 = \sum_{k=0}^n \binom{n+1}{k} B_k
+    .. math :: n+1 = \sum_{k=0}^n \binom{n+1}{k} B_k
 
     They are also commonly defined by their exponential generating
-    function, which is `\frac{x}{e^x - 1}`. For odd indices > 1, the
-    Bernoulli numbers are zero.
+    function, which is `\frac{x}{1 - e^{-x}}`. For odd indices > 1,
+    the Bernoulli numbers are zero.
 
     The Bernoulli polynomials satisfy the analogous formula:
 
-    .. math :: B_n(x) = \sum_{k=0}^n \binom{n}{k} B_k x^{n-k}
+    .. math :: B_n(x) = \sum_{k=0}^n (-1)^k \binom{n}{k} B_k x^{n-k}
 
     Bernoulli numbers and Bernoulli polynomials are related as
-    `B_n(0) = B_n`.
+    `B_n(1) = B_n`.
+
+    The generalized Bernoulli function `\operatorname{B}(s, a)`
+    is defined for any complex `s` and `a`, except where `a` is a
+    nonpositive integer and `s` is not a nonnegative integer. It is
+    an entire function of `s` for fixed `a`, related to the Hurwitz
+    zeta function by
+
+    .. math:: \operatorname{B}(s, a) = \begin{cases}
+              -s \zeta(1-s, a) & s \ne 0 \\ 1 & s = 0 \end{cases}
+
+    When `s` is a nonnegative integer this function reduces to the
+    Bernoulli polynomials: `\operatorname{B}(n, x) = B_n(x)`. When
+    `a` is omitted it is assumed to be 1, yielding the (ordinary)
+    Bernoulli function which interpolates the Bernoulli numbers and is
+    related to the Riemann zeta function.
 
     We compute Bernoulli numbers using Ramanujan's formula:
 
@@ -452,21 +467,36 @@ class bernoulli(Function):
     .. math :: S(n) = \sum_{k=1}^{[n/6]} \binom{n+3}{n-6k} B_{n-6k}
 
     This formula is similar to the sum given in the definition, but
-    cuts 2/3 of the terms. For Bernoulli polynomials, we use the
-    formula in the definition.
+    cuts `\frac{2}{3}` of the terms. For Bernoulli polynomials, we use
+    the formula in the definition.
+
+    For `n` a nonnegative integer and `s`, `a`, `x` arbitrary complex numbers,
 
     * ``bernoulli(n)`` gives the nth Bernoulli number, `B_n`
+    * ``bernoulli(s)`` gives the Bernoulli function `\operatorname{B}(s)`
     * ``bernoulli(n, x)`` gives the nth Bernoulli polynomial in `x`, `B_n(x)`
+    * ``bernoulli(s, a)`` gives the generalized Bernoulli function
+      `\operatorname{B}(s, a)`
+
+    .. versionchanged:: 1.12
+        ``bernoulli(1)`` gives `+\frac{1}{2}` instead of `-\frac{1}{2}`.
+        This choice of value confers several theoretical advantages [5]_,
+        including the extension to complex parameters described above
+        which this function now implements. The previous behavior, defined
+        only for nonnegative integers `n`, can be obtained with
+        ``(-1)**n*bernoulli(n)``.
 
     Examples
     ========
 
     >>> from sympy import bernoulli
-
+    >>> from sympy.abc import x
     >>> [bernoulli(n) for n in range(11)]
-    [1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66]
+    [1, 1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66]
     >>> bernoulli(1000001)
     0
+    >>> bernoulli(3, x)
+    x**3 - 3*x**2/2 + x/2
 
     See Also
     ========
@@ -480,6 +510,10 @@ class bernoulli(Function):
     .. [2] https://en.wikipedia.org/wiki/Bernoulli_polynomial
     .. [3] http://mathworld.wolfram.com/BernoulliNumber.html
     .. [4] http://mathworld.wolfram.com/BernoulliPolynomial.html
+    .. [5] Peter Luschny, "The Bernoulli Manifesto",
+           http://luschny.de/math/zeta/The-Bernoulli-Manifesto.html
+    .. [6] Peter Luschny, "An introduction to the Bernoulli function",
+           https://arxiv.org/abs/2009.06743
 
     """
 
@@ -513,12 +547,11 @@ def eval(cls, n, x=None):
         elif n.is_zero:
             return S.One
         elif n.is_integer is False or n.is_nonnegative is False:
-            raise ValueError("Bernoulli numbers are defined only"
-                             " for nonnegative integer indices.")
+            return
         # Bernoulli numbers
         elif x is None:
             if n is S.One:
-                return -S.Half
+                return S.Half
             elif n.is_odd and (n-1).is_positive:
                 return S.Zero
             elif n.is_Number:
@@ -1269,7 +1302,7 @@ class genocchi(Function):
 
     >>> from sympy import genocchi, Symbol
     >>> [genocchi(n) for n in range(1, 9)]
-    [1, -1, 0, 1, 0, -3, 0, 17]
+    [-1, -1, 0, 1, 0, -3, 0, 17]
     >>> n = Symbol('n', integer=True, positive=True)
     >>> genocchi(2*n + 1)
     0
diff --git a/sympy/functions/special/zeta_functions.py b/sympy/functions/special/zeta_functions.py
index de9a39b37e..71fdb6a34e 100644
--- a/sympy/functions/special/zeta_functions.py
+++ b/sympy/functions/special/zeta_functions.py
@@ -510,9 +510,7 @@ def eval(cls, s, a=None):
         sint = s.is_Integer
         if a is None:
             a = S.One
-        if s.is_zero:
-            return S.Half - a
-        elif sint and s.is_nonpositive:
+        if sint and s.is_nonpositive:
             return bernoulli(1-s, a) / (s-1)
         elif a is S.One:
             if sint and s.is_even:
diff --git a/sympy/printing/numpy.py b/sympy/printing/numpy.py
index 5ab29585a1..5f6186e411 100644
--- a/sympy/printing/numpy.py
+++ b/sympy/printing/numpy.py
@@ -400,6 +400,10 @@ def _print_airybiprime(self, expr):
                 self._module_format("scipy.special.airy"),
                 self._print(expr.args[0]))
 
+    def _print_bernoulli(self, expr):
+        # scipy's bernoulli is inconsistent with SymPy's so rewrite
+        return self._print(expr._eval_rewrite_as_zeta(*expr.args))
+
     def _print_Integral(self, e):
         integration_vars, limits = _unpack_integral_limits(e)
 


Test patch:
diff --git a/sympy/functions/combinatorial/tests/test_comb_numbers.py b/sympy/functions/combinatorial/tests/test_comb_numbers.py
index 4c726b081b..98e1265781 100644
--- a/sympy/functions/combinatorial/tests/test_comb_numbers.py
+++ b/sympy/functions/combinatorial/tests/test_comb_numbers.py
@@ -41,7 +41,7 @@ def test_carmichael():
 
 def test_bernoulli():
     assert bernoulli(0) == 1
-    assert bernoulli(1) == Rational(-1, 2)
+    assert bernoulli(1) == Rational(1, 2)
     assert bernoulli(2) == Rational(1, 6)
     assert bernoulli(3) == 0
     assert bernoulli(4) == Rational(-1, 30)
@@ -72,7 +72,30 @@ def test_bernoulli():
     assert isinstance(bernoulli(2 * l + 1), bernoulli)
     assert isinstance(bernoulli(2 * m + 1), bernoulli)
     assert bernoulli(2 * n + 1) == 0
-    raises(ValueError, lambda: bernoulli(-2))
+
+    assert bernoulli(x, 1) == bernoulli(x)
+
+    assert str(bernoulli(0.0, 2.3).evalf(n=10)) == '1.000000000'
+    assert str(bernoulli(1.0).evalf(n=10)) == '0.5000000000'
+    assert str(bernoulli(1.2).evalf(n=10)) == '0.4195995367'
+    assert str(bernoulli(1.2, 0.8).evalf(n=10)) == '0.2144830348'
+    assert str(bernoulli(1.2, -0.8).evalf(n=10)) == '-1.158865646 - 0.6745558744*I'
+    assert str(bernoulli(3.0, 1j).evalf(n=10)) == '1.5 - 0.5*I'
+    assert str(bernoulli(I).evalf(n=10)) == '0.9268485643 - 0.5821580598*I'
+    assert str(bernoulli(I, I).evalf(n=10)) == '0.1267792071 + 0.01947413152*I'
+    assert bernoulli(x).evalf() == bernoulli(x)
+
+
+def test_bernoulli_rewrite():
+    from sympy.functions.elementary.piecewise import Piecewise
+    n = Symbol('n', integer=True, nonnegative=True)
+
+    assert bernoulli(-1).rewrite(zeta) == pi**2/6
+    assert bernoulli(-2).rewrite(zeta) == 2*zeta(3)
+    assert not bernoulli(n, -3).rewrite(zeta).has(harmonic)
+    assert bernoulli(-4, x).rewrite(zeta) == 4*zeta(5, x)
+    assert isinstance(bernoulli(n, x).rewrite(zeta), Piecewise)
+    assert bernoulli(n+1, x).rewrite(zeta) == -(n+1) * zeta(-n, x)
 
 
 def test_fibonacci():
@@ -425,7 +448,7 @@ def test_catalan():
 
 
 def test_genocchi():
-    genocchis = [1, -1, 0, 1, 0, -3, 0, 17]
+    genocchis = [-1, -1, 0, 1, 0, -3, 0, 17]
     for n, g in enumerate(genocchis):
         assert genocchi(n + 1) == g
 
diff --git a/sympy/utilities/tests/test_lambdify.py b/sympy/utilities/tests/test_lambdify.py
index caa80eac84..f06c210c01 100644
--- a/sympy/utilities/tests/test_lambdify.py
+++ b/sympy/utilities/tests/test_lambdify.py
@@ -11,6 +11,7 @@
 from sympy.core.singleton import S
 from sympy.core.symbol import (Dummy, symbols)
 from sympy.functions.combinatorial.factorials import (RisingFactorial, factorial)
+from sympy.functions.combinatorial.numbers import bernoulli
 from sympy.functions.elementary.complexes import Abs
 from sympy.functions.elementary.exponential import exp, log
 from sympy.functions.elementary.hyperbolic import acosh
@@ -1467,6 +1468,14 @@ def test_scipy_special_math():
     assert abs(cm1(1e-20) + 5e-41) < 1e-200
 
 
+def test_scipy_bernoulli():
+    if not scipy:
+        skip("scipy not installed")
+
+    bern = lambdify((x,), bernoulli(x), modules='scipy')
+    assert bern(1) == 0.5
+
+
 def test_cupy_array_arg():
     if not cupy:
         skip("CuPy not installed")