diff --git a/acb.h b/acb.h
index 4c88a9d77..db9f27834 100644
--- a/acb.h
+++ b/acb.h
@@ -780,7 +780,9 @@ void acb_lgamma(acb_t y, const acb_t x, slong prec);
 void acb_log_sin_pi(acb_t res, const acb_t z, slong prec);
 void acb_digamma(acb_t y, const acb_t x, slong prec);
 void acb_zeta(acb_t z, const acb_t s, slong prec);
+void acb_bernoulli(acb_t z, const acb_t s, slong prec);
 void acb_hurwitz_zeta(acb_t z, const acb_t s, const acb_t a, slong prec);
+void acb_bernoulli_gen(acb_t z, const acb_t s, const acb_t a, slong prec);
 void acb_polygamma(acb_t res, const acb_t s, const acb_t z, slong prec);
 
 void acb_bernoulli_poly_ui(acb_t res, ulong n, const acb_t x, slong prec);
diff --git a/acb/zeta.c b/acb/zeta.c
index 61f435024..a874d2007 100644
--- a/acb/zeta.c
+++ b/acb/zeta.c
@@ -43,3 +43,14 @@ acb_zeta(acb_t z, const acb_t s, slong prec)
     acb_dirichlet_zeta(z, s, prec);
 }
 
+void
+acb_bernoulli_gen(acb_t z, const acb_t s, const acb_t a, slong prec)
+{
+    acb_dirichlet_bernoulli_gen(z, s, a, prec);
+}
+
+void
+acb_bernoulli(acb_t z, const acb_t s, slong prec)
+{
+    acb_dirichlet_bernoulli(z, s, prec);
+}
diff --git a/acb_dirichlet.h b/acb_dirichlet.h
index d79eb2772..f056b1478 100644
--- a/acb_dirichlet.h
+++ b/acb_dirichlet.h
@@ -42,11 +42,13 @@ void acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K);
 void acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec);
 void acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec);
 void acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec);
+void acb_dirichlet_bernoulli(acb_t res, const acb_t s, slong prec);
 
 void acb_dirichlet_zeta_jet_rs(acb_ptr res, const acb_t s, slong len, slong prec);
 void acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, slong prec);
 
 void acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec);
+void acb_dirichlet_bernoulli_gen(acb_t res, const acb_t s, const acb_t a, slong prec);
 
 void acb_dirichlet_lerch_phi_integral(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec);
 void acb_dirichlet_lerch_phi_direct(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec);
diff --git a/acb_dirichlet/hurwitz.c b/acb_dirichlet/hurwitz.c
index 8f6d890d7..b3c7a13a3 100644
--- a/acb_dirichlet/hurwitz.c
+++ b/acb_dirichlet/hurwitz.c
@@ -55,3 +55,35 @@ acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec)
     _acb_poly_zeta_cpx_series(res, s, a, 0, 1, prec);
 }
 
+void
+acb_dirichlet_bernoulli_gen(acb_t res, const acb_t s, const acb_t a, slong prec)
+{
+    if (acb_is_zero(s))
+    {
+        acb_one(res);
+        return;
+    }
+
+    if (acb_is_int(s) && arf_sgn(arb_midref(acb_realref(s))) > 0 &&
+        arf_cmpabs_ui(arb_midref(acb_realref(s)), prec / 2) < 0)
+    {
+        slong n = arf_get_si(arb_midref(acb_realref(s)), ARF_RND_FLOOR);
+
+        acb_bernoulli_poly_ui(res, n, a, prec);
+        return;
+    }
+    else
+    {
+        acb_t t;
+
+        acb_init(t);
+
+        acb_neg(t, s);
+        acb_add_ui(res, t, 1, prec);
+        acb_dirichlet_hurwitz(res, res, a, prec);
+        acb_mul(res, t, res, prec);
+
+        acb_clear(t);
+    }
+}
+
diff --git a/acb_dirichlet/zeta.c b/acb_dirichlet/zeta.c
index 6c671e3f5..cb5f725aa 100644
--- a/acb_dirichlet/zeta.c
+++ b/acb_dirichlet/zeta.c
@@ -92,3 +92,26 @@ acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec)
     }
 }
 
+void
+acb_dirichlet_bernoulli(acb_t res, const acb_t s, slong prec)
+{
+    if (acb_is_zero(s))
+    {
+        acb_one(res);
+        return;
+    }
+    else
+    {
+        acb_t t;
+
+        acb_init(t);
+
+        acb_neg(t, s);
+        acb_add_ui(res, t, 1, prec);
+        acb_dirichlet_zeta(res, res, prec);
+        acb_mul(res, t, res, prec);
+
+        acb_clear(t);
+    }
+}
+
diff --git a/arb.h b/arb.h
index 2e40f576d..9ef38795e 100644
--- a/arb.h
+++ b/arb.h
@@ -567,7 +567,9 @@ void arb_gamma_fmpq(arb_t z, const fmpq_t x, slong prec);
 void arb_gamma_fmpz(arb_t z, const fmpz_t x, slong prec);
 void arb_digamma(arb_t y, const arb_t x, slong prec);
 void arb_zeta(arb_t z, const arb_t s, slong prec);
+void arb_bernoulli(arb_t z, const arb_t s, slong prec);
 void arb_hurwitz_zeta(arb_t z, const arb_t s, const arb_t a, slong prec);
+void arb_bernoulli_gen(arb_t z, const arb_t s, const arb_t a, slong prec);
 
 void arb_rising_ui(arb_t z, const arb_t x, ulong n, slong prec);
 void arb_rising_fmpq_ui(arb_t y, const fmpq_t x, ulong n, slong prec);
diff --git a/arb/hurwitz_zeta.c b/arb/hurwitz_zeta.c
index 9b7916037..6ac0f9d7c 100644
--- a/arb/hurwitz_zeta.c
+++ b/arb/hurwitz_zeta.c
@@ -40,3 +40,26 @@ arb_hurwitz_zeta(arb_t res, const arb_t s, const arb_t z, slong prec)
     }
 }
 
+void
+arb_bernoulli_gen(arb_t res, const arb_t s, const arb_t z, slong prec)
+{
+    if (arb_is_zero(s))
+    {
+        arb_one(res);
+        return;
+    }
+    else
+    {
+        arb_t t;
+
+        arb_init(t);
+
+        arb_neg(t, s);
+        arb_add_ui(res, t, 1, prec);
+        arb_hurwitz_zeta(res, res, z, prec);
+        arb_mul(res, t, res, prec);
+
+        arb_clear(t);
+    }
+}
+
diff --git a/arb/zeta.c b/arb/zeta.c
index ac865ef49..2f29c92c5 100644
--- a/arb/zeta.c
+++ b/arb/zeta.c
@@ -23,3 +23,26 @@ arb_zeta(arb_t y, const arb_t s, slong prec)
     acb_clear(t);
 }
 
+void
+arb_bernoulli(arb_t res, const arb_t s, slong prec)
+{
+    if (arb_is_zero(s))
+    {
+        arb_one(res);
+        return;
+    }
+    else
+    {
+        arb_t t;
+
+        arb_init(t);
+
+        arb_neg(t, s);
+        arb_add_ui(res, t, 1, prec);
+        arb_zeta(res, res, prec);
+        arb_mul(res, t, res, prec);
+
+        arb_clear(t);
+    }
+}
+
diff --git a/doc/source/acb.rst b/doc/source/acb.rst
index 0aaf50633..1805d6d3e 100644
--- a/doc/source/acb.rst
+++ b/doc/source/acb.rst
@@ -1015,6 +1015,18 @@ Zeta function
 
     This is a wrapper of :func:`acb_dirichlet_hurwitz`.
 
+.. function:: void acb_bernoulli(acb_t z, const acb_t s, slong prec)
+
+    Sets *z* to the value of the Bernoulli function `B(s)`.
+
+    This is a wrapper of :func:`acb_dirichlet_bernoulli`.
+
+.. function:: void acb_bernoulli_gen(acb_t z, const acb_t s, const acb_t a, slong prec)
+
+    Sets *z* to the value of the generalized Bernoulli function `B(s, a)`.
+
+    This is a wrapper of :func:`acb_dirichlet_bernoulli_gen`.
+
 .. function:: void acb_bernoulli_poly_ui(acb_t res, ulong n, const acb_t x, slong prec)
 
     Sets *res* to the value of the Bernoulli polynomial `B_n(x)`.
diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst
index ecd5c60c5..81b1b02f0 100644
--- a/doc/source/acb_dirichlet.rst
+++ b/doc/source/acb_dirichlet.rst
@@ -132,6 +132,11 @@ Riemann zeta function
     `\xi(s) = \frac{1}{2} s (s-1) \pi^{-s/2} \Gamma(\frac{1}{2} s) \zeta(s)`.
     The functional equation for xi is `\xi(1-s) = \xi(s)`.
 
+.. function:: void acb_dirichlet_bernoulli(acb_t res, const acb_t s, slong prec)
+
+    Sets *res* to the Bernoulli function `B(s) = -s \zeta(1-s)`, with
+    the limiting value of 1 used when `s = 0` [Lus2020]_.
+
 Riemann-Siegel formula
 -------------------------------------------------------------------------------
 
@@ -207,6 +212,11 @@ Hurwitz zeta function
     when `a = 1`. Some other special cases may also be handled by direct
     formulas. In general, Euler-Maclaurin summation is used.
 
+.. function:: void acb_dirichlet_bernoulli_gen(acb_t res, const acb_t s, const acb_t a, slong prec)
+
+    Computes the generalized Bernoulli function `B(s, a) = -s \zeta(1-s, a)`, with
+    the limiting value of 1 used when `s = 0` [Lus2020]_.
+
 Hurwitz zeta function precomputation
 -------------------------------------------------------------------------------
 
diff --git a/doc/source/arb.rst b/doc/source/arb.rst
index 6a60f7c79..3591553dc 100644
--- a/doc/source/arb.rst
+++ b/doc/source/arb.rst
@@ -1455,6 +1455,11 @@ Zeta function
     For computing derivatives with respect to `s`,
     use :func:`arb_poly_zeta_series`.
 
+.. function:: void arb_bernoulli(arb_t z, const arb_t s, slong prec)
+
+    Sets *z* to the value of the Bernoulli function `B(s) = -s \zeta(1-s)`,
+    with the limiting value of 1 used when `s = 0` [Lus2020]_.
+
 .. function:: void arb_hurwitz_zeta(arb_t z, const arb_t s, const arb_t a, slong prec)
 
     Sets *z* to the value of the Hurwitz zeta function `\zeta(s,a)`.
@@ -1462,6 +1467,12 @@ Zeta function
     For computing derivatives with respect to `s`,
     use :func:`arb_poly_zeta_series`.
 
+.. function:: void arb_bernoulli_gen(arb_t z, const arb_t s, const arb_t a, slong prec)
+
+    Sets *z* to the value of the generalized Bernoulli function
+    `B(s,a) = -s \zeta(1-s,a)`, with the limiting value of 1 used when
+    `s = 0` [Lus2020]_.
+
 Bernoulli numbers and polynomials
 -------------------------------------------------------------------------------
 
diff --git a/doc/source/credits.rst b/doc/source/credits.rst
index 1636e98e6..5a1662788 100644
--- a/doc/source/credits.rst
+++ b/doc/source/credits.rst
@@ -249,6 +249,8 @@ Bibliography
 
 .. [Leh1970] \R. S. Lehman, "On the Distribution of Zeros of the Riemann Zeta-Function", Proc. of the London Mathematical Society 20(3) (1970), 303-320, https://doi.org/10.1112/plms/s3-20.2.303
 
+.. [Lus2020] \P. Luschny, "An introduction to the Bernoulli function", preprint (2020), https://arxiv.org/abs/2009.06743
+
 .. [Mic2007] \N. Michel, "Precise Coulomb wave functions for a wide range of complex l, eta and z", Computer Physics Communications, Volume 176, Issue 3, (2007), 232-249, https://doi.org/10.1016/j.cpc.2006.10.004
 
 .. [Miy2010] \S. Miyajima, "Fast enclosure for all eigenvalues in generalized eigenvalue problems", Journal of Computational and Applied Mathematics, 233 (2010), 2994-3004, https://dx.doi.org/10.1016/j.cam.2009.11.048