add arb_bernoulli_fmpz for completeness

arb.h

@@ -717,6 +717,7 @@ void arb_zeta_ui_vec_odd(arb_ptr x, ulong start, slong num, slong prec);
 void arb_zeta_ui_vec(arb_ptr x, ulong start, slong num, slong prec);
 void arb_bernoulli_ui(arb_t b, ulong n, slong prec);
 void arb_bernoulli_ui_zeta(arb_t b, ulong n, slong prec);
+void arb_bernoulli_fmpz(arb_t b, const fmpz_t n, slong prec);
 
 void arb_bernoulli_poly_ui(arb_t res, ulong n, const arb_t x, slong prec);
 

arb/bernoulli_fmpz.c

@@ -0,0 +1,73 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb.h"
+
+void
+arb_bernoulli_fmpz(arb_t res, const fmpz_t n, slong prec)
+{
+    if (fmpz_cmp_ui(n, UWORD_MAX) <= 0)
+    {
+        if (fmpz_sgn(n) >= 0)
+            arb_bernoulli_ui(res, fmpz_get_ui(n), prec);
+        else
+            arb_zero(res);
+    }
+    else if (fmpz_is_odd(n))
+    {
+        arb_zero(res);
+    }
+    else
+    {
+        arb_t t;
+        slong wp;
+
+        arb_init(t);
+        wp = prec + 2 * fmpz_bits(n);
+
+        /* zeta(n) ~= 1 */
+        arf_one(arb_midref(res));
+        mag_one(arb_radref(res));
+        mag_mul_2exp_si(arb_radref(res), arb_radref(res), WORD_MIN);
+
+        /* |B_n| = 2 * Gamma(n)! / (2*pi)^n * zeta(n) */
+        arb_gamma_fmpz(t, n, wp);
+        arb_mul_fmpz(t, t, n, wp);
+        arb_mul(res, res, t, wp);
+
+        arb_const_pi(t, wp);
+        arb_mul_2exp_si(t, t, 1);
+        arb_pow_fmpz(t, t, n, wp);
+
+        arb_div(res, res, t, prec);
+        arb_mul_2exp_si(res, res, 1);
+
+        if (fmpz_fdiv_ui(n, 4) == 0)
+            arb_neg(res, res);
+
+        arb_clear(t);
+    }
+}
+

doc/source/arb.rst

@@ -1218,11 +1218,21 @@ Bernoulli numbers and polynomials
 
 .. function:: void arb_bernoulli_ui(arb_t b, ulong n, slong prec)
 
-    Sets `b` to the numerical value of the Bernoulli number `B_n` accurate
-    to *prec* bits, computed by a division of the exact fraction if `B_n` is in
-    the global cache or the exact numerator roughly is larger than
-    *prec* bits, and using :func:`arb_bernoulli_ui_zeta`
-    otherwise. This function reads `B_n` from the global cache
+.. function:: void arb_bernoulli_fmpz(arb_t b, const fmpz_t n, slong prec)
+
+    Sets `b` to the numerical value of the Bernoulli number `B_n`
+    approximated to *prec* bits.
+
+    The internal precision is increased automatically to give an accurate
+    result. Note that, with huge *fmpz* input, the output will have a huge
+    exponent and evaluation will accordingly be slower.
+
+    A single division from the exact fraction of `B_n` is used if this value
+    is in the global cache or the exact numerator roughly is larger than
+    *prec* bits. Otherwise, the Riemann zeta function is used
+    (see :func:`arb_bernoulli_ui_zeta`).
+
+    This function reads `B_n` from the global cache
     if the number is already cached, but does not automatically extend
     the cache by itself.
 
@@ -1233,9 +1243,8 @@ Bernoulli numbers and polynomials
 
     To avoid potential infinite recursion, we explicitly call the
     Euler product implementation of the zeta function.
-    We therefore assume that the precision is small
-    enough and `n` large enough for the Euler product to converge
-    rapidly (otherwise this function will effectively hang).
+    This method will only give high accuracy if the precision is small
+    enough compared to `n` for the Euler product to converge rapidly.
 
 .. function:: void arb_bernoulli_poly_ui(arb_t res, ulong n, const arb_t x, slong prec)