Generalized zeta function imprecise for large negative s

[Expander]

The current implementation of the generalized zeta function $\zeta(s,z)$ suffers from an imprecision for large negative $s$:

julia> SpecialFunctions.zeta(-100, 2.0)
-1.921642939869165e44

The correct result should be $-1$, because for negative integer $s$ one would have
$$\zeta(-n,a)=-\frac{B_{n+1}(a)}{n+1}$$
[Wikipedia]
For example Mathematica gives the correct result:

In[]:= Zeta[-100, 2] // N[#,17]&                                            
Out[]= -1.0000000000000000
In[]:= -BernoulliB[100 + 1, 2]/(100 + 1)
Out[]= -1