Entire (generalised) Bernoulli function

[Parcly-Taxel]

This PR is like sympy/sympy#23984, but for Arb. It implements the generalised Bernoulli function from @PeterLuschny's "An introduction to the Bernoulli function":

$$B(s,a)=\begin{cases}1&s=0\\-s\zeta(1-s,a)&s\ne0\end{cases}$$

as well as the ordinary version $B(s)=B(s,1)$. These functions, like the Riemann/Hurwitz zeta functions, come in arb and acb versions too.

Note that the integer-only Bernoulli number functions (arb_bernoulli_ui, etc.) are not touched, and that this PR is inapplicable to FLINT since it has no Hurwitz zeta function.

It might be possible to improve accuracy near the zeta pole (i.e. the origin, $s=0$) by directly using the Maclaurin series expansion of the Bernoulli function (which is entire for fixed $a$), but this is not implemented here.

---

Example of usage:

#include "arb.h"
#include "acb.h"
#include <stdio.h>

int main()
{
    arb_t s, a, x;
    arb_init(s);
    arb_init(a);
    arb_init(x);
    arb_set_d(s, 0.9);
    arb_set_d(a, 0.7);
    arb_bernoulli(x, s, 100);
    arb_printn(x, 50, 0);
    printf("\n");
    arb_bernoulli_gen(x, s, a, 100);
    arb_printn(x, 50, 0);
    printf("\n");
    arb_clear(s);
    arb_clear(a);
    arb_clear(x);

    acb_t s2, a2, x2;
    acb_init(s2);
    acb_init(a2);
    acb_init(x2);
    arb_set_d(acb_realref(s2), 0.4);
    arb_set_d(acb_imagref(s2), -1.2);
    arb_set_d(acb_realref(a2), -0.7);
    arb_set_d(acb_imagref(a2), 0.3);
    acb_bernoulli(x2, s2, 100);
    acb_printn(x2, 50, 0);
    printf("\n");
    acb_bernoulli_gen(x2, s2, a2, 100);
    acb_printn(x2, 50, 0);
    printf("\n");
    acb_clear(s2);
    acb_clear(a2);
    acb_clear(x2);
}

Output:

[0.54273376787061753404885841995 +/- 4.12e-30]
[0.248587684454566258768404567803 +/- 5.42e-31]
[0.66727675711496067968005244886 +/- 2.95e-30] + [0.627508782028046109043402027631 +/- 8.85e-31]*I
[30.583734710361609686614823596 +/- 2.22e-28] + [25.559309022539927788601797491 +/- 5.14e-28]*I