It is tough to say for sure. The formula above is a rough approximation of the number of summands needed to calculate the function. I think simply changing the "15" in that formula to use config.epsilon would work just fine. The formula written mathematically is:
To calculate $\zeta(a+ib)$ to $d$ decimals places of accuracy, set $n=1.3d+0.9|b|$. I don't see a major issue with just setting $d$ to the big number precision.

Hmm, I gave it a try but it doesn't magically result in the desired precision (also not with Complex numbers).

You may have better luck using an alternative error bound. I have calculated $$n=\frac{\ln(3\cdot 10^d(2|t|+1)s^{|t|\pi/2})}{\ln(3+\sqrt{2})}$$. If this is still troublesome we may have to resort to the full error bound involving the gamma function which isn't too friendly :(

The error bound for the formulation is
$$\displaystyle|\gamma_n(s)| \leq \frac{2}{(3+\sqrt{8})^n} \frac{1}{|\Gamma(s)|}\frac{1}{|1-2^{1-s}|} \leq \frac{3}{3+\sqrt{8}} \frac{(1+2|t|)e^{|t|\pi / 2}}{|1-2^{1-s}|}$$
and so what I did was use the last expression and take the limit as the real part of $s$ goes to infinity to get the bound needed for large real values. The function is asymptotic anyway so it works fine. I have tested this using both the limit variation and standard one using $s$ and both produce pratically the exact same value for n with Re(s) that are not very close to 1. If they are close to 1 we can just make use of the full formula for calculating n.
$$\displaystyle n=\frac{\ln\left(3\cdot\frac{10^{d}\left(2\left|t\right|+1\right)e^{\frac{t\pi}{2}}}{|1-2^{1-(\sigma + it)}|}\right)}{\ln\left(3+\sqrt{2}\right)}$$
and for values far larger than 1 just use the limit and say that the term $|1-2^{1-(\sigma + it)}|$ goes to 1
After crunching the numbers, I think you could use
$$0.74+1.55d\ +\ 0.67\ln\left(2|t|+1\right)+1.06|t|$$
The original $1.3d + 0.9|t|$ was actually an underestimate of the actual error. This new one is an overestimate which can be good or bad depending on how you look at it. The full error (if you choose to use it) would be this in expanded form
$$n=\frac{1}{\ln\left(3+\sqrt{2}\right)}\left(\ln\left(3\right)+d\ln\left(10\right)+\ln\left(2\left|t\right|+1\right)+\frac{\left|t\right|\pi}{2}-\ln\left(\left|1-2^{\left(1-\left(a+it\right)\right)}\right|\right)\right)$$
Let me know if this results in the desired precision or not and I can help troubleshoot

Thanks for working out these formulas for the bounds. I did a bit of trial and error but no success yet. I tried the following for the determineDigits argument for numbers for example (and the same but with Math.abs(s.im) for complex numbers):

value => 0.74 + 1.55 * (12 + 1) + 0.67 * Math.log(2 * Math.abs(value) + 1) + 1.06 * Math.abs(value)

Help troubleshooting this would be very welcome.

This is very strange indeed. When I tested it in python it seemed to work fine with all the precision. I suspect an issue with the d(k,n) function possibly when n gets large the factorials become too large to store precisely causing loss in decimal digits. Maybe replacing it with a big number variant may work...?

Maybe replacing it with a big number variant may work...?

I don't think that will make a difference since the number of digits is not a huge number nor does it need a lot of precision.

Any other ideas to get the precision of the function under control?

@gwhitney do you maybe have an idea on how to ensure the accuracy of this new zeta function?

If not, I propose we leave it as it is and add a comment in the docs that the function has a limited precision of about 6 digits.

Sorry for the long delay, I was on vacation. I'm not sure why the precision is not working but I'll continue searching. For now we can just leave it at 6 digits of accuracy I think

I've tried to extract the internal zetaNumeric function so we can use it for the number only mathjs, but it is not turning out how I want. I'll leave it as it is for now. We can always have a look again in the future.

Published now in v11.10.0