Derive zeta values of even integers from the Euler-Maclaurin formula. Euler showed:
\begin{equation}
 B_{2 k} = (-1)^{k+1} \frac{2 \,  (2 \, k)!}{ (2 \, \pi)^{2 k}} 
   \zeta(2 k)
\end{equation}
for $k=1,2, \cdots$. We could from here find $\zeta(2k)$ in terms of the
even Bernoulli coefficients $B_{2k}$.
How can we derive the equivalent representation by using the
Euler Maclaurin formula ?
Thanks.
Update
Here is what I have done:
\begin{eqnarray*}
  \zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s} \quad , \quad  \mathrm{Re}(s) > 1.
\end{eqnarray*}
Euler used the Euler-Maclaurin series to find values of the Riemann Zeta function.
  We want to use  the 
Euler-Maclaurin formula with $f(x)=1/x^s=x^{-s}$. We know that 
\begin{eqnarray*}
  f^{(0)}(x) &=& \frac{1}{x^s} \\
  f^{(1)}(x) &=& -\frac{s}{x^{s+1}} \\
  f^{(2)}(x) &=& \frac{s(s+1)}{x^{s+2}} \\
  &\vdots& \\
  f^{(i)}(x) &=& (-1)^{i+1} \frac{s(s+1) \dots (s+i-2)}{x^{s+i}}  =
  (-1)^{i+1}   \frac{\Gamma(s+i)}{ \Gamma(s)} \frac{1}{ x^{s+i}}.
\end{eqnarray*}
We then write using $h=1$, $a=1$, $b=\infty$
\begin{eqnarray*}
  \sum_{i=1}^{\infty} \frac{1}{n^s} = \int_1^{\infty} \frac{dx}{x^s} + \frac{1}{2}  +
  \left . \sum_{i=1}^{m} \frac{B_{2i}}{(2i)!} 
  \frac{\Gamma(s+2 i-1)}{ \Gamma(s)} \frac{1}{ x^{s+2 i-1}} \right |_1^{\infty} + R_{2m}
\end{eqnarray*}
with 
\begin{eqnarray*}
  R_{2m} =
  -\int_1^{\infty} \mathrm{B}_{2m} 
  \left \{ x-1 \right \} 
  \frac{\Gamma(s+2 m)}{(2m)! \, \Gamma(s)} \frac{dx}{x^{s+2m}}.
\end{eqnarray*}
That is
\begin{eqnarray}
  \zeta(s) =  \frac{1}{s-1} + \frac{1}{2} -
   \sum_{i=1}^{m} \frac{B_{2i}}{(2i)!} 
  \frac{\Gamma(s+2 i-1)}{ \Gamma(s)}  + R_{2m}
  \label{zetazeta}
\end{eqnarray}
This is an important equation since it establishes an analytic continuation for the 
$\zeta(s)$ function. We observe that the first fraction is an analytic function except
for $s=1$, then the sum of quotients of $\Gamma$  functions is analytic except for
isolated singularities in the negative integer arguments. Finally, since
the Bernoulli polynomial $B_{2m}\{x-1\}$ is bounded the residual is a convergent integral for $s + 2m >1$, so we can extend the 
convergence as far as $s > 1-2m$, for any positive number $m$. 
The question here is what "$m$" to choose. If I choose $m=1$, I did the computations and found something which made no sense. 
 A: Should you accept to use the functional equation then let's start with the  Euler Maclaurin formula for zeta in the specific case where the binomial sum is convergent (see for example this thread) :
$$\tag{1}\zeta(s)= \sum_{k=1}^N \frac 1{k^s} {+\frac 1{(s-1)\;N^{s-1}}} + \sum_{i=1}^{\infty} \frac{B_{i}\;s^{(i-1)}}{i!\;N^{s-1+i}},\quad N\in\mathbb{N}^+$$ 
with $s^{(i)}$ the "rising factorial" : $\;s^{(i)}=s(s+1)\cdots(s+k-1),\,s^{(0)}=1$ and $B_1=-\dfrac 12$.
(I used the general formula $(1)$ but your formula supposing $N=1$ should give the same result!)
For $\;s:=1-j\;$ with $\,j\,$ a positive integer the EMl formula is applied to the polynomial $\,f(k):=k^{1-j}\;$ so that $(1)$ is valid will a finite sum at the right from $\;(1-j)^{(j)}=0\;$ giving : 
\begin{align}
\zeta(1-j)&= \sum_{k=1}^N {k^{j-1}} -\frac {N^{j}}{j\;} + \sum_{i=1}^{j} \frac{B_{i}\;(1-j)^{(i-1)}N^{j-i}}{i!\;}\\
\zeta(1-j)&= \sum_{k=1}^N {k^{j-1}} -\frac {N^{j}}{j\;} + \sum_{i=1}^{j} (-1)^{i-1}\frac{(j-1)!}{i!\,(j-i)!}B_{i}\;N^{j-i}\\
\tag{2}\zeta(1-j)&= \sum_{k=1}^N {k^{j-1}} -\frac {N^{j}}{j\;} - \frac 1j\sum_{i=1}^{j} (-1)^{i}{j\choose{i}}B_{i}\;N^{j-i}\\
\end{align}
But the Faulhaber formula gives us :
$$\tag{3}\sum_{k=1}^N {k^{j-1}}= {1\over{j}}\sum_{i=0}^{j-1} (-1)^i {j\choose{i}} B_i\; N^{j-i}$$
When we combine these two expressions all the terms of the sum disappear except for $i=0$ and $i=j\,$ and :
$$\zeta(1-j)= {1\over{j}}(-1)^0 {j\choose{0}} B_0\; N^{j} -\frac {N^{j}}{j\;} - \frac 1j(-1)^{j}{j\choose{j}}B_{j}\;N^{0}= -(-1)^{j}\frac {B_{j}}j$$
Since $\,B_{2n+1}=0\,$ for any integer $n>1$ this becomes the simple :
$$\tag{4}\zeta(1-j)= -\frac {B_{j}}j$$
To conclude we will use the promised functional equation applied to $\,j:=1-2m$ :
\begin{align}
\zeta(1-2m) &= 2^{1-2m}\pi^{-2m}\ \sin\left(\frac{\pi (1-2m)}2\right)\ \Gamma(2m)\ \zeta(2m)\\
-\frac {B_{2m}}{2m}&= 2(2\pi)^{-2m}\ \cos\left(\pi m\right)\ \Gamma(2m)\;\zeta(2m) \\
B_{2m}&= 2(2\pi)^{-2m}\ (-1)^{m+1}\ \Gamma(2m+1)\;\zeta(2m) \\
\end{align}
with the wished result :
$$\tag{5}\boxed{\displaystyle B_{2m}= (-1)^{m+1}\,\frac{2\,(2m)!}{(2\pi)^{2m}}\;\zeta(2m) }$$ 
$$-$$
References to Euler's derivation of the Euler–Maclaurin formula are here.
