Questions about iterating the Euler-Maclaurin summation formula

Question

Introduction

The Euler–Maclaurin summation formula is as follows for a positive integer $p$ and a continuous function $f(\cdot)$ that is $p$ times continuously differentiable on the interval $[m,n]$ : $$\sum_{i=m}^{n} f(i) = \int_{m}^{n} f(x) \, dx + \frac{f(n)+f(m)}{2} + \sum_{k=1}^{\lfloor p/2 \rfloor} \frac{B_{2k}}{(2k)!}\big{(}f^{(2k-1)}(n)-f^{(2k-1)}(m)\big{)} + R_{p}.$$

Here, the remainder term satisfies the inequality $$R_{p} \leq \frac{2\zeta(p)}{(2 \pi)^{p}} \int_{m}^{n} \big{|}f^{(p)} (x) \big{|} \, dx .$$

Usually, this formula is used to approximate the value of sums by picking some suitable finite value of $p$. However, I'm interested in ascertaining whether it could be used to find closed forms of different classes of sums. To that end, let $p \to \infty$ and assume that $$ \lim_{p \to \infty} R_{p} = 0 $$ for the function $f$ we're looking into.

Now, we are only left with three terms. The third term is the one that often prevents the sum from being evaluated exactly. The idea I have in mind is to apply the Euler-Maclaurin summation formula again to this term.

Iterating the E-M formula: an attempt at an example

Suppose $f(x) = 1/x^{3}$, $m=1$ and $n \to \infty$. The first and second terms of the E-M formula are easy to evaluate in this case, so let's focus on the third one. There is no obvious way to find this sum. So what I try is to iterate over the E-M formula by setting $$g_{a}(x) = \frac{B_{2x}}{(2x)!} f^{(2x-1)}(a) $$ and try to find \begin{align} \lim_{n \to \infty} \Bigg{(} \sum_{k=1}^{\infty} \frac{B_{2k}}{(2k)!}\big{(}f^{(2k-1)}(n)-f^{(2k-1)}(1)\big{)} \Bigg{)} &= \int_{1}^{\infty} \big{(}g_{\infty}(x) - g_{1}(x)\big{)} \, dx \\ & \quad + \Bigg{(} \frac{g_{\infty}(x)-g_{1}(x)}{2} \Bigg{)} \Bigg{|}^{x=\infty}_{x=1} \\ & \quad + \sum_{m=1}^{\infty} \frac{B_{2m}}{(2m)!} \Big{(} g^{(2m-1)}_{\infty}(x) - g^{(2m-1)}_{1}(x) \Big{)} \bigg{|}^{x=\infty}_{x=1} \qquad \quad (*) \end{align}

My aim to obtain all three terms on the right side of $(*)$. Below, I focus on the first term: the definite integral.

Now, we know that $$B_{2x} = \frac{2 \Gamma(2x+1)}{(2\pi)^{2x}} \zeta(2x) .$$ Moreover, we have $$(2x)! = \Gamma(2x+1) .$$ Finally, we derive the last part by means of the Riemann-Liouville fractional derivative: \begin{align} \frac{d^{2x-1}}{dy^{2x-1}} y^{-3} &= \frac{d^{2x}}{dy^{2x}} \Big{(} - \frac{1}{2} y^{-2} \Big{)} \\ &= \Big{(}-\frac{1}{2}\Big{)} \cdot 6 \frac{d^{x}}{dy^{x}} y^{-4} \\ &= \Big{(}-\frac{1}{2}\Big{)} \cdot 6 \cdot (-1)^{x} \frac{\Gamma(4+x)}{\Gamma(4)} y^{-(4+x)} \\ &= \Big{(}- \frac{1}{2}e^{i \pi x} \Big{)} \Gamma(4+x) y^{-(4+x)}. \end{align} Combining these expressions, we obtain $$g_{a}(x) = - \frac{e^{i \pi x}\Gamma(4+x)a^{-(4+x)}}{(2 \pi)^{2x}} \zeta(2x). $$

As $a \to \infty$, this expression amounts to zero (because $x$ is positive). But we also have $$g_{1}(x) = - \frac{e^{i \pi x}\Gamma(4+x)}{(2 \pi)^{2x}} \zeta(2x). $$

If we wish to iterate the E-M formula and find, for instance, the first term, we must integrate this expression from $x=1$ to $x=\infty$. This is where I get worried, because the term $\Gamma(4+x)$ grows much quicker than $(2 \pi)^{2x}$, so I'm afraid this integral doesn't converge.

Questions

1. Is my derivation of $g_{a}(x)$ and its corresponding E-M formula correct? I am particularly interested in whether I did the Riemann-Liouville fractional derivative part without errors. If not, could you show what it should be?
2. Has an iterated approach of the Euler-Maclaurin formula been formulated already in earlier work?
3. If the Riemann-Liouville fractional derivative doesn't work for this application, is there another type of fractional derivative that does allow for iterating the E-M formula? (Considering there are many different definitions of the fractional derivative of a function).

Note: earlier, I asked a slightly less elaborate version of this question on MSE.