s np
>>> from scipy.special import bernoulli, zeta
>>> bernoulli(4)
array([ 1.        , -0.5       ,  0.16666667,  0.        , -0.03333333])

The Wikipedia article ([2]_) points out the relationship between the
Bernoulli numbers and the zeta function, ``B_n^+ = -n * zeta(1 - n)``
for ``n > 0``:

>>> n = np.arange(1, 5)
>>> -n * zeta(1 - n)
array([ 0.5       ,  0.16666667, -0.        , -0.03333333])

Note that, in the notation used in the wikipedia article,
`bernoulli` computes ``B_n^-`` (i.e. it used the convention that
``B_1`` is -1/2).  The relation given above is for ``B_n^+``, so the
sign of 0.5 does not match the output of ``bernoulli(4)``.

r