Repository: mpmath/mpmath
Pull request #724 (2023-10-07 08:10:31)

Add plus flag to select the B_1 sign convention for bernoulli/bernfrac
Based on #639 // @Parcly-Taxel

Some links:
- https://github.com/sympy/sympy/issues/23866 (the root of)
- https://groups.google.com/g/sage-devel/c/G4sbKoibXK4/m/nruZdm9xBwAJ
- https://groups.google.com/g/sage-devel/c/G4sbKoibXK4/m/LmO04RclCAAJ (Fredrik opinion)
- https://github.com/sagemath/sage/issues/34521
- https://github.com/flintlib/flint2/pull/1179
- https://github.com/gnu-octave/symbolic/issues/1217

I'm not sure if we should deprecate `-0.5` convention.  Does it help, for example, to SymPy (which already uses `+0.5` choice)?  Not so much, I think.



Patch:
diff --git a/docs/functions/numtheory.rst b/docs/functions/numtheory.rst
index e7d51af..35222c8 100644
--- a/docs/functions/numtheory.rst
+++ b/docs/functions/numtheory.rst
@@ -18,11 +18,11 @@ Bernoulli numbers and polynomials
 
 :func:`~mpmath.bernoulli`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
-.. autofunction:: mpmath.bernoulli(n)
+.. autofunction:: mpmath.bernoulli
 
 :func:`~mpmath.bernfrac`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-.. autofunction:: mpmath.bernfrac(n)
+.. autofunction:: mpmath.bernfrac
 
 :func:`~mpmath.bernpoly`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/mpmath/ctx_fp.py b/mpmath/ctx_fp.py
index b93d0e8..7494bc0 100644
--- a/mpmath/ctx_fp.py
+++ b/mpmath/ctx_fp.py
@@ -1,4 +1,5 @@
 import cmath
+import functools
 import math
 import sys
 
@@ -18,7 +19,6 @@ def __init__(ctx):
 
         # Override SpecialFunctions implementation
         ctx.loggamma = libfp.loggamma
-        ctx._bernoulli_cache = {}
         ctx.pretty = False
 
         ctx._init_aliases()
@@ -56,12 +56,9 @@ def f_wrapped(ctx, *args, **kwargs):
         f_wrapped.__doc__ = function_docs.__dict__.get(name, f.__doc__)
         setattr(cls, name, f_wrapped)
 
-    def bernoulli(ctx, n):
-        cache = ctx._bernoulli_cache
-        if n in cache:
-            return cache[n]
-        cache[n] = to_float(mpf_bernoulli(n, 53, 'n'), strict=True)
-        return cache[n]
+    @functools.lru_cache
+    def bernoulli(ctx, n, plus=False):
+        return to_float(mpf_bernoulli(n, 53, 'n', plus=plus), strict=True)
 
     pi = libfp.pi
     e = math.e
diff --git a/mpmath/ctx_mp.py b/mpmath/ctx_mp.py
index 34a0942..8334c0d 100644
--- a/mpmath/ctx_mp.py
+++ b/mpmath/ctx_mp.py
@@ -219,8 +219,8 @@ def _agm(ctx, a, b=1):
         else: b = b._mpc_
         return ctx.make_mpc(libmp.mpc_agm(a, b, prec, rounding))
 
-    def bernoulli(ctx, n):
-        return ctx.make_mpf(libmp.mpf_bernoulli(int(n), *ctx._prec_rounding))
+    def bernoulli(ctx, n, plus=False):
+        return ctx.make_mpf(libmp.mpf_bernoulli(int(n), *ctx._prec_rounding, plus=plus))
 
     def _zeta_int(ctx, n):
         return ctx.make_mpf(libmp.mpf_zeta_int(int(n), *ctx._prec_rounding))
diff --git a/mpmath/function_docs.py b/mpmath/function_docs.py
index 0c9dc7c..1a3cad3 100644
--- a/mpmath/function_docs.py
+++ b/mpmath/function_docs.py
@@ -2228,6 +2228,9 @@
 returns a floating-point approximation. To obtain an exact
 fraction, use :func:`~mpmath.bernfrac` instead.
 
+Optional ``plus`` flag (default: False) control the sign choice of
+the `B_1` value (default: `-0.5`).
+
 **Examples**
 
 Numerical values of the first few Bernoulli numbers::
@@ -2282,6 +2285,11 @@
 
 For larger `n`, `B_n` is evaluated in terms of the Riemann zeta
 function.
+
+**References**
+
+1. https://en.wikipedia.org/wiki/Bernoulli_number
+
 """
 
 stieltjes = r"""
diff --git a/mpmath/libmp/gammazeta.py b/mpmath/libmp/gammazeta.py
index 9c086ff..8e3ca4a 100644
--- a/mpmath/libmp/gammazeta.py
+++ b/mpmath/libmp/gammazeta.py
@@ -380,7 +380,7 @@ def bernoulli_size(n):
 
 BERNOULLI_PREC_CUTOFF = bernoulli_size(MAX_BERNOULLI_CACHE)
 
-def mpf_bernoulli(n, prec, rnd=None):
+def mpf_bernoulli(n, prec, rnd=None, plus=False):
     """Computation of Bernoulli numbers (numerically)"""
     if n < 2:
         if n < 0:
@@ -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 if plus else mpf_neg(fhalf)
     # For odd n > 1, the Bernoulli numbers are zero
     if n & 1:
         return fzero
@@ -466,13 +466,16 @@ def mpf_bernoulli_huge(n, prec, rnd=None):
         v = mpf_neg(v)
     return mpf_pos(v, prec, rnd or round_fast)
 
-def bernfrac(n):
+def bernfrac(n, plus=False):
     r"""
     Returns a tuple of integers `(p, q)` such that `p/q = B_n` exactly,
     where `B_n` denotes the `n`-th Bernoulli number. The fraction is
     always reduced to lowest terms. Note that for `n > 1` and `n` odd,
     `B_n = 0`, and `(0, 1)` is returned.
 
+    Optional ``plus`` flag (default: False) control the sign choice of
+    the `B_1` value (default: `(-1, 2)`).
+
     **Examples**
 
     The first few Bernoulli numbers are exactly::
@@ -540,10 +543,12 @@ def bernfrac(n):
     2. The Bernoulli Number Page:
        http://www.bernoulli.org/
 
+    3. https://en.wikipedia.org/wiki/Bernoulli_number
+
     """
     n = int(n)
     if n < 3:
-        return [(1, 1), (-1, 2), (1, 6)][n]
+        return [(1, 1), (1 if plus else -1, 2), (1, 6)][n]
     if n & 1:
         return (0, 1)
     q = 1


Test patch:
diff --git a/mpmath/tests/test_fp.py b/mpmath/tests/test_fp.py
index 104fa9e..b0a6a5d 100644
--- a/mpmath/tests/test_fp.py
+++ b/mpmath/tests/test_fp.py
@@ -21,6 +21,8 @@
 
 """
 
+import pytest
+
 from mpmath import fp
 
 
@@ -122,12 +124,13 @@ def test_fp_expj():
     assert ae(fp.expjpi(0.75), (-0.7071067811865475244 + 0.7071067811865475244j))
     assert ae(fp.expjpi(2+3j), (0.000080699517570304599239 + 0.0j))
 
-def test_fp_bernoulli():
-    assert ae(fp.bernoulli(0), 1.0)
-    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)
+@pytest.mark.parametrize('plus', [True, False])
+def test_fp_bernoulli(plus):
+    assert ae(fp.bernoulli(0, plus), 1.0)
+    assert ae(fp.bernoulli(1, plus), 0.5 if plus else -0.5)
+    assert ae(fp.bernoulli(2, plus), 0.16666666666666666667)
+    assert ae(fp.bernoulli(10, plus), 0.075757575757575757576)
+    assert ae(fp.bernoulli(11, plus), 0.0)
 
 def test_fp_gamma():
     assert ae(fp.gamma(1), 1.0)
diff --git a/mpmath/tests/test_gammazeta.py b/mpmath/tests/test_gammazeta.py
index dd0c352..bdeeb31 100644
--- a/mpmath/tests/test_gammazeta.py
+++ b/mpmath/tests/test_gammazeta.py
@@ -1,3 +1,5 @@
+import pytest
+
 from mpmath import (altzeta, apery, barnesg, bell, bernfrac, bernoulli,
                     bernpoly, beta, binomial, catalan, digamma, e, euler,
                     eulerpoly, fac, fac2, factorial, fadd, ff, findroot, fp,
@@ -11,48 +13,49 @@
 def test_zeta_int_bug():
     assert mpf_zeta_int(0, 10) == from_float(-0.5)
 
-def test_bernoulli():
-    assert bernfrac(0) == (1,1)
-    assert bernfrac(1) == (-1,2)
-    assert bernfrac(2) == (1,6)
-    assert bernfrac(3) == (0,1)
-    assert bernfrac(4) == (-1,30)
-    assert bernfrac(5) == (0,1)
-    assert bernfrac(6) == (1,42)
-    assert bernfrac(8) == (-1,30)
-    assert bernfrac(10) == (5,66)
-    assert bernfrac(12) == (-691,2730)
-    assert bernfrac(18) == (43867,798)
-    p, q = bernfrac(228)
+@pytest.mark.parametrize('plus', [True, False])
+def test_bernoulli(plus):
+    assert bernfrac(0, plus) == (1,1)
+    assert bernfrac(1, plus) == (1,2) if plus else (-1,2)
+    assert bernfrac(2, plus) == (1,6)
+    assert bernfrac(3, plus) == (0,1)
+    assert bernfrac(4, plus) == (-1,30)
+    assert bernfrac(5, plus) == (0,1)
+    assert bernfrac(6, plus) == (1,42)
+    assert bernfrac(8, plus) == (-1,30)
+    assert bernfrac(10, plus) == (5,66)
+    assert bernfrac(12, plus) == (-691,2730)
+    assert bernfrac(18, plus) == (43867,798)
+    p, q = bernfrac(228, plus)
     assert p % 10**10 == 164918161
     assert q == 625170
-    p, q = bernfrac(1000)
+    p, q = bernfrac(1000, plus)
     assert p % 10**10 == 7950421099
     assert q == 342999030
     mp.dps = 15
-    assert bernoulli(0) == 1
-    assert bernoulli(1) == -0.5
-    assert bernoulli(2).ae(1./6)
-    assert bernoulli(3) == 0
-    assert bernoulli(4).ae(-1./30)
-    assert bernoulli(5) == 0
-    assert bernoulli(6).ae(1./42)
-    assert str(bernoulli(10)) == '0.0757575757575758'
-    assert str(bernoulli(234)) == '7.62772793964344e+267'
-    assert str(bernoulli(10**5)) == '-5.82229431461335e+376755'
-    assert str(bernoulli(10**8+2)) == '1.19570355039953e+676752584'
+    assert bernoulli(0, plus) == 1
+    assert bernoulli(1, plus) == 0.5 if plus else -0.5
+    assert bernoulli(2, plus).ae(1./6)
+    assert bernoulli(3, plus) == 0
+    assert bernoulli(4, plus).ae(-1./30)
+    assert bernoulli(5, plus) == 0
+    assert bernoulli(6, plus).ae(1./42)
+    assert str(bernoulli(10, plus)) == '0.0757575757575758'
+    assert str(bernoulli(234, plus)) == '7.62772793964344e+267'
+    assert str(bernoulli(10**5, plus)) == '-5.82229431461335e+376755'
+    assert str(bernoulli(10**8+2, plus)) == '1.19570355039953e+676752584'
 
     mp.dps = 50
-    assert str(bernoulli(10)) == '0.075757575757575757575757575757575757575757575757576'
-    assert str(bernoulli(234)) == '7.6277279396434392486994969020496121553385863373331e+267'
-    assert str(bernoulli(10**5)) == '-5.8222943146133508236497045360612887555320691004308e+376755'
-    assert str(bernoulli(10**8+2)) == '1.1957035503995297272263047884604346914602088317782e+676752584'
+    assert str(bernoulli(10, plus)) == '0.075757575757575757575757575757575757575757575757576'
+    assert str(bernoulli(234, plus)) == '7.6277279396434392486994969020496121553385863373331e+267'
+    assert str(bernoulli(10**5, plus)) == '-5.8222943146133508236497045360612887555320691004308e+376755'
+    assert str(bernoulli(10**8+2, plus)) == '1.1957035503995297272263047884604346914602088317782e+676752584'
 
     mp.dps = 1000
-    assert bernoulli(10).ae(mpf(5)/66)
+    assert bernoulli(10, plus).ae(mpf(5)/66)
 
     mp.dps = 50000
-    assert bernoulli(10).ae(mpf(5)/66)
+    assert bernoulli(10, plus).ae(mpf(5)/66)
 
     mp.dps = 15