I wanted to know if there's any existing algorithm to compute exact values of the Riemann Zeta function at even natural numbers? For example, it should compute $\zeta(4)$ as exactly $\frac{\pi^4}{90}$ i.e.in fraction form, not a decimal approximation.
If there is such an algorithm, could you state its complexity as well? Similarly, is there an algorithm to calculate the values of the Dirichlet Beta function at odd natural numbers?
On
The value is known to be: $$ \zeta(2 k) = \frac{(-1)^{k + 1} 4^k \pi^{2 k} B_{2 k}}{2 (2 k)!} $$ The Bernoulli numbers page probably has the best algorithms to compute the $B_{2 k}$.
In addition to @mixedmath's comment: just invert the matrix $P_n$ of the size $nxn$, for example $n=6$ $$ P_n =P_6 = \small \left[ \begin{array} {} 1 & . & . & . & . & . \\ 1 & 2 & . & . & . & . \\ 1 & 3 & 3 & . & . & . \\ 1 & 4 & 6 & 4 & . & . \\ 1 & 5 & 10 & 10 & 5 & . \\ 1 & 6 & 15 & 20 & 15 & 6 \end{array} \right] $$ which is simply a submatrix of the Patscalmatrix reduced by the diagonal. You get $$ G_6=P_6^{-1} =\small \begin{bmatrix} 1 & . & . & . & . & . \\ -1/2 & 1/2 & . & . & . & . \\ 1/6 & -1/2 & 1/3 & . & . & . \\ 0 & 1/4 & -1/2 & 1/4 & . & . \\ -1/30 & 0 & 1/3 & -1/2 & 1/5 & . \\ 0 & -1/12 & 0 & 5/12 & -1/2 & 1/6 \end{bmatrix} $$ and in the first column you'll find the Bernoulli-numbers. You can do this for any size $n$. Becasue the matrix is triangular, the inversion is really a simple task.
After you have the Bernoulli-numbers, you can express the zetas at even indexes by the well known conversion formula (due to L. Euler) as rational multiple of powers of $\pi$