Numerical evaluation of Hurwitz zeta function
Note that your variables names are uncommon, I will use the standard names from DLMF $$\zeta(s,a) = \sum\limits_{n=0}^{\infty}\frac{1}{(n+a)^s}$$ A good numerical algorithm uses the Euler-Maclaurin summation (DLMF or Wolfram) in the form $$\zeta(s,a) = \sum\limits_{k=0}^{n}\frac{1}{(a+k)^s} + \frac{(a+n)^{1-s}}{s-1} + \sum\limits_{k=0}^{\infty} \frac{B_{k+1}}{(k+1)!}\: \frac{(s+k-1)_k}{(a+n)^{s+k}}$$ where $(\cdot)_k$ is the Pochhammer symbol. Using the properties of the Bernoulli numbers $B_1=-1/2\;$ and $B_{2k+1}=0\;$ this can be rewritten as $$\zeta(s,a) = \sum\limits_{k=0}^{n}\frac{1}{(a+k)^s} + \frac{(a+n)^{1-s}}{s-1} - \frac{1}{2(a+n)^{s}} + \sum\limits_{k=1}^{\infty} \frac{B_{2k}}{(2k)!}\: \frac{(s+2k-2)_{2k-1}}{(a+n)^{s+2k-1}} \cdot $$ The upper bound of the first sum should be selected depending on your arithmetic (e.g $n=8$ or $n=9$ for double etc). The second (infinite) sum can computed by summing the terms until there is no change (i.e. the current term is smaller in magnitude than current sum).
This algorithm is efficient for the primary range $s>0, a>0:\;$ For 80-bit floating point arithmetic I get $\zeta(2.5,0.75)\approx 2.49154238551193522\;$ with $10$ terms of the second sum, and $\zeta(10,0.25)\approx 1048576.10768311475\;$ you only need two terms.
For negative $s$ you can use Bernoulli polynomials (for $s=-1,-2\dots$) or the reflection formula from DLMF.
Edit: The number of terms depends mainly on $s$ and $n$. Here some intermediate values from my implementation for 53-bit IEEE double. For the given $(s,a)\;$ pairs $n=9.\;$ s1 is the first sum from $k=0\dots n,\;$ s2 = s1 + the next two single terms, and s2 the total sum after finishing the last sum.
ZetaH(10, 0.25)
s1 = 1048576.10768311471
s2 = 1048576.10768311494
s3 = 1048576.10768311494 1 term
ZetaH(5.0, 0.25)
s1 = 1024.34894710160870
s2 = 1024.34897386673833
s3 = 1024.34897452658106 6 terms
ZetaH(2.5, 0.75)
s1 = 2.47125718157242114
s2 = 2.49147059743736321
s3 = 2.49154238551193474 8 terms
ZetaH(0.75, 0.25)
s1 = 6.13122575577889695
s2 = -0.938864666031481665
s3 = -0.937591964187151961 9 terms
In actual implementations, you can adjust the $n$ value, and of course there are other speed-up methods: Terminate the first sum, if there are no changes, pre-compute a certain numbers of the complicated coefficients of the second sum etc.
You find a basic C implementation in the file zeta.c of the well known Cephes library http://www.moshier.net/double.zip
Edit 2: Here the data for the two sums. The first column shows n (the upper limit of the first sum), j is the index of the second sum, for which the convergence criterion is achieved. The second sum is also terminated, if the terms are non-decreasing, this is the case for n=4,5 and the last three s,a pairs, in the table it is shown as ---.
s,a= s,a= s,a= s,a=
10,0.25 5,0.25 2.5,0.75 0.75,0.25
n=4 j=7 --- --- ---
n=5 j=3 --- --- ---
n=6 j=1 j=8 j=11 j=12
n=7 j=1 j=6 j=9 j=10
n=8 j=0 j=5 j=8 j=8
n=9 j=0 j=5 j=7 j=8
n=10 j=4 j=7 j=7
n=11 j=4 j=6 j=7
n=12 j=4 j=6 j=6
n=13 j=3 j=6 j=6
n=14 j=3 j=5 j=6
n=15 j=3 j=5 j=6
n=16 j=3 j=5 j=5
n=17 j=3 j=5 j=5
n=18 j=3 j=5 j=5
n=19 j=2 j=4 j=5
n=20 j=2 j=4 j=5
n=21 j=2 j=4 j=5
n=22 j=2 j=4 j=5
n=23 j=2 j=4 j=5
n=24 j=2 j=4 j=5
n=25 j=2 j=4 j=4
n=26 j=2 j=4 j=4
The first formula given by @gammatester agrees with that given on the Wolfram site, but the second, rewritten, form seems to be in error.
The first, simple, error is that the exponent for the $(a+n)$ term in the denominator of Bernoulli series is wrong by 2 - should be $s+2k-1$, but a more involved error is in interpreting the Pochhammer symbol when summing over the even indices, $2k$. As given, in the rewritten formula it amounts to $(s)2k+1$. This is incorrect: consider first the Pochhammer form in the original form as $(s + k-1)k$. This means the first factor in the Pochhammer symbol advances by one each time, and the last one by 2. In the even sum, it should therefore be $(s + 2k-2)2k-1$. The first term is $(s)1 =s$, the second is $(s+2)(s+3)(s+4)$, etc. The last factor is $s+4k-4$.
Numerically the quickest way to evaluate in the loop is to increment $k$ by one, two times. $(s+k-1)k$, is advanced one step by dividing by the first factor then adding two more.
I have implemented the correct algorithm in Fortran and Mathematica (to higher than double precision), and am in agreement with @allard that the error term after the first explicit summation - i.e. with the series including Bernoulli numbers - does not decrease as expected. The error saturates long before the series terms begin to diverge, and is not a precision issue.