A derivation of the Euler-Maclaurin formula?

Question

The generating function for the Bernoulli numbers $B_n$ is $$\frac{x}{e^x-1}=\sum_{n=0}^\infty\frac{B_n}{n!}x^n$$

The sum of an infinite geometric series is $$\frac{1}{1-x}=\sum_{k=0}^\infty x^k$$

Replacing $x$ with $e^x$ yields $$\frac{1}{1-e^x}=\sum_{k=0}^\infty(e^x)^k$$

Multiplying both sides by $-x$ yields $$\frac{x}{e^x-1}=-\sum_{k=0}^\infty(e^x)^kx$$

Hence $$\sum_{k=0}^\infty(e^x)^kx+\sum_{n=0}^\infty\frac{B_n}{n!}x^n=0$$

Replacing $x$ with the derivative operator $D$ yields $$\sum_{k=0}^\infty(e^D)^kD+\sum_{n=0}^\infty\frac{B_n}{n!}D^n=0$$

The exponential of the derivative operator is the shift operator $T$: $$e^D=T$$

Hence $$\sum_{k=0}^\infty T^kD+\sum_{n=0}^\infty\frac{B_n}{n!}D^n=0$$

Applying this to a function $f$ yields $$\sum_{k=0}^\infty (T^kDf)(x)+\sum_{n=0}^\infty\frac{B_n}{n!}(D^nf)(x)=0$$

which is equivalent by $(Tf)(x)=f(x+1)$ to $$\sum_{k=0}^\infty (Df)(x+k)+\sum_{n=0}^\infty\frac{B_n}{n!}(D^nf)(x)=0$$

If $Df=g$, this means $$\sum_{k=0}^\infty g(x+k)+\int_a^xg(t)\,\mathrm{d}t+\sum_{n=0}^\infty\frac{B_{n+1}}{(n+1)!}(D^ng)(x)=0$$

Since the values of the Riemann zeta function at the negative integers are $$\zeta(-n)=-\frac{B_{n+1}}{(n+1)}$$

This can be expressed as $$\sum_{k=0}^\infty g(x+k)+\int_a^xg(t)\,\mathrm{d}t-\sum_{n=0}^\infty\frac{\zeta(-n)}{n!}(D^ng)(x)=0$$

or $$\sum_{k=0}^\infty g(x+k)=\sum_{n=0}^\infty\frac{\zeta(-n)}{n!}(D^ng)(x)-\int_a^xg(t)\,\mathrm{d}t$$

Have I derived a version of the Euler-Maclaurin formula correctly? Is an additional error term needed somewhere? Why does the expression on the right side of the equation seem to depend on $a$, while the one on the left does not?

How can I manipulate this procedure to obtain the more traditional version of the formula, e.g. with finite sums? I know, for example, that one can express a finite sum of powers using

$$\frac{1-x^{n+1}}{1-x}=1+x+x^2+x^3+\ldots+x^n=\sum_{k=0}^n x^k$$