Chowla and Hartung provide an "algorithm" for computing Bernoulli numbers in this paper.
In particular, if the Bernoulli numbers are defined by $$\frac{x}{e^x-1}=1-\frac{x}{2}+\sum_{n=1}^\infty \frac{(-1)^{n-1}B_n x^{2n}}{(2n)!},$$ then for $m\ge 1$: $$2(2^{2m}-1)B_m=\lfloor \varphi_m\rfloor +1$$ where $$\varphi_m=\frac{2(2^{2m}-1)(2m)!}{2^{2m-1}\pi^{2m}}\sum_{n=1}^{3m}\frac{1}{n^{2m}}.$$
There is one problem though. When trying to implement this "algorithm", I realized that I need some sequence of rational approximations to $\pi$. With increasing $m$, I need better and better rational approximations to $\pi$. This is because, for large $m$, $\pi^{2m}$ and $(\pi\pm\varepsilon)^{2m}$ are very different even for small $\varepsilon$.
The question is: To what accuracy do I need to approximate $\pi$ in order for this "algorithm" to work for a given $m$?
This problem is not discussed in the paper at all so I posted this question.
This question has also been asked on MSE (with no answers so far).