I was playing around with ways to calculate Bernoulli numbers; for this post I will take their generating function as $x/(1-e^{-x})$, that is, $$ \sum_{n=0}^{\infty} \frac{B_n}{n!}x^n = \frac{x}{1-e^{-x}} $$Here I want $B_1=\textbf{+}1/2$ but it shouldn't really matter. Computing $B_n$ this way often involves L'Hopital's Rule, which I was hoping to avoid.
One way around this is to use the Maclaurin polynomial for $e^{-x}$: to calculate $B_n$, one could replace $e^{-x}$ with its degree $n+1$ Maclaurin polynomial, differentiate $n$ times, and evaluate at $x=0$. For example, to calculate $B_4$, we would put $$ B_4=\left.\frac{d^4}{dx^4} \frac{x}{1-\left(1-x+x^2/2-x^3/6+x^4/24-x^5/120\right)}\right|_{x=0} $$ $$ =\left.(14400 \left(7 x^{12}-105 x^{11}+826 x^{10}-4130 x^9+12755 x^8-21700 x^7-4700 x^6+129720 x^5-324800 x^4+432000 x^3-192000 x^2-144000 x-57600\right)\left(x^4-5 x^3+20 x^2-60 x+120\right)^{-5}\right|_{x=0} $$ $$ =\frac{-14400\cdot 57600}{120^5}=\frac{-1}{30} $$If we're game for derivatives of rational functions, we should have faith that this procedure should work since $e^{-x}$ is analytic in a neighborhood of zero.
What I tried to do next doesn't immediately seem like it should work and is the reason for my question. Instead of expanding $e^{-x}$ as a Taylor series, what if we expanded $1/(1-e^{-x})$ as a geometric series? $$ \frac{x}{1-e^{-x}} = x \sum_{k=0}^{\infty} e^{-kx}; $$one immediate problem is that equality only holds for $\Re(x)>0$. Nevertheless, I decided to truncate this series at $r$ terms and differentiate $n$ times. The derivatives termwise are easy: here are the results for $n=1,\ldots,6$: $$ \left\{r+1,-r (r+1),\frac{1}{2} r (r+1) (2 r+1),-r^2 (r+1)^2,\frac{1}{6} r (r+1) (2 r+1) \left(3 r^2+3 r-1\right),-\frac{1}{2} r^2 (r+1)^2 \left(2 r^2+2 r-1\right)\right\} $$Now if we integrate these over $-1\le r\le 0$, somehow we get the Bernoulli numbers back: $$ \left\{\frac{1}{2},\frac{1}{6},0,-\frac{1}{30},0,\frac{1}{42}\right\} $$My question: why is this justified and why does this work?