Change bernoulli(1) to +1/2

This commit is part of the project at
https://github.com/sympy/sympy/pull/23926

Co-authored-by: Sergey B Kirpichev <skirpichev@gmail.com>

Adapted from #639

mpmath/function_docs.py

@@ -2228,6 +2228,10 @@ The Bernoulli numbers are rational numbers, but this function
 returns a floating-point approximation. To obtain an exact
 fraction, use :func:`~mpmath.bernfrac` instead.
 
+Note that `B_1=+0.5`; this choice of value confers several theoretical
+advantages [1]. The previous behavior can be obtained with
+``(-1)**n*bernoulli(n)``.
+
 **Examples**
 
 Numerical values of the first few Bernoulli numbers::
@@ -2238,7 +2242,7 @@ Numerical values of the first few Bernoulli numbers::
     ...     print("%s %s" % (n, bernoulli(n)))
     ...
     0 1.0
-    1 -0.5
+    1 0.5
     2 0.166666666666667
     3 0.0
     4 -0.0333333333333333
@@ -2282,6 +2286,12 @@ guaranteed to be fast.
 
 For larger `n`, `B_n` is evaluated in terms of the Riemann zeta
 function.
+
+**References**
+
+1. P. Luschny, "The Bernoulli Manifesto",
+   https://luschny.de/math/zeta/The-Bernoulli-Manifesto.html
+
 """
 
 stieltjes = r"""

mpmath/libmp/gammazeta.py

@@ -388,7 +388,7 @@ def mpf_bernoulli(n, prec, rnd=None):
         if n == 0:
             return fone
         if n == 1:
-            return mpf_neg(fhalf)
+            return fhalf
     # For odd n > 1, the Bernoulli numbers are zero
     if n & 1:
         return fzero
@@ -473,6 +473,10 @@ def bernfrac(n):
     always reduced to lowest terms. Note that for `n > 1` and `n` odd,
     `B_n = 0`, and `(0, 1)` is returned.
 
+    Note that `B_1=+\frac{1}{2}`; this choice of value confers several
+    theoretical advantages [3]. The previous behavior can be obtained by
+    multiplying ``p`` by ``(-1)**n``.
+
     **Examples**
 
     The first few Bernoulli numbers are exactly::
@@ -483,7 +487,7 @@ def bernfrac(n):
         ...     print("%s %s/%s" % (n, p, q))
         ...
         0 1/1
-        1 -1/2
+        1 1/2
         2 1/6
         3 0/1
         4 -1/30
@@ -540,10 +544,13 @@ def bernfrac(n):
     2. The Bernoulli Number Page:
        http://www.bernoulli.org/
 
+    3. P. Luschny, "The Bernoulli Manifesto",
+       https://luschny.de/math/zeta/The-Bernoulli-Manifesto.html
+
     """
     n = int(n)
     if n < 3:
-        return [(1, 1), (-1, 2), (1, 6)][n]
+        return [(1, 1), (1, 2), (1, 6)][n]
     if n & 1:
         return (0, 1)
     q = 1

mpmath/tests/test_fp.py

@@ -124,7 +124,7 @@ def test_fp_expj():
 
 def test_fp_bernoulli():
     assert ae(fp.bernoulli(0), 1.0)
-    assert ae(fp.bernoulli(1), -0.5)
+    assert ae(fp.bernoulli(1), 0.5)
     assert ae(fp.bernoulli(2), 0.16666666666666666667)
     assert ae(fp.bernoulli(10), 0.075757575757575757576)
     assert ae(fp.bernoulli(11), 0.0)

mpmath/tests/test_gammazeta.py

@@ -13,7 +13,7 @@ def test_zeta_int_bug():
 
 def test_bernoulli():
     assert bernfrac(0) == (1,1)
-    assert bernfrac(1) == (-1,2)
+    assert bernfrac(1) == (1,2)
     assert bernfrac(2) == (1,6)
     assert bernfrac(3) == (0,1)
     assert bernfrac(4) == (-1,30)
@@ -31,7 +31,7 @@ def test_bernoulli():
     assert q == 342999030
     mp.dps = 15
     assert bernoulli(0) == 1
-    assert bernoulli(1) == -0.5
+    assert bernoulli(1) == 0.5
     assert bernoulli(2).ae(1./6)
     assert bernoulli(3) == 0
     assert bernoulli(4).ae(-1./30)