ag.algebraic geometry - Euler-Maclaurin formula and Riemann-Roch

Let $Df$ denote the derivative of a function $f(x)$ and $bigtriangledown f=f(x)-f(x-1)$ be the discrete derivative. Using the Taylor series expansion for $f(x-1)$, we easily get $bigtriangledown = 1- e^{-D}$ or, by taking the inverses,
$$ frac{1}{bigtriangledown} = frac{1}{1-e^{-D}} =
frac{1}{D}cdot frac{D}{1-e^{-D}}=
frac{1}{D} + frac12+
sum_{k=1}^{infty} B_{2k}frac{D^{2k-1}}{(2k)!}
,$$
where $B_{2k}$ are Bernoulli numbers.

(Edit: I corrected the signs to adhere to the most common conventions.)

Here, $(1/D)g$ is the opposite to the derivative, i.e. the integral; adding the limits this becomes a definite integral $int_0^n g(x)dx$. And $(1/bigtriangledown)g$ is the opposite to the discrete derivative, i.e. the sum $sum_{x=1}^n g(x)$. So the above formula, known as Euler-Maclaurin formula, allows one, sometimes, to compute the discrete sum by using the definite integral and some error terms.

Usually, there is a nontrivial remainder in this formula. For example, for $g(x)=1/x$, the remainder is Euler's constant $gammasimeq 0.57$. Estimating the remainder and analyzing the convergence of the power series is a long story, which is explained for example in the nice book "Concrete Mathematics" by Graham-Knuth-Patashnik. But the power series becomes finite with zero remainder if $g(x)$ is a polynomial. OK, so far I am just reminding elementary combinatorics.

Now, for my question. In the (Hirzebruch/Grothendieck)-Riemann-Roch formula one of the main ingredients is the Todd class which is defined as the product, going over Chern roots $alpha$, of the expression $frac{alpha}{1-e^{-alpha}}$. This looks so similar to the above, and so suggestive (especially because in the Hirzebruch's version
$$chi(X,F) = h^0(F)-h^1(F)+dots = int_X ch(F) Td(T_X)$$
there is also an "integral", at least in the notation) that it makes me wonder: is there a connection?

The obvious case to try (which I did) is the case when $X=mathbb P^n$ and $F=mathcal O(d)$. But the usual proof in that case is a residue computation which, to my eye, does not look anything like Euler-Maclaurin formula.

But is there really a connection?

---

An edit after many answers: Although the connection with Khovanskii-Pukhlikov's paper and the consequent work, pointed out by Dmitri and others, is undeniable, it is still not obvious to me how the usual Riemann-Roch for $X=mathbb P^n$ and $F=mathcal O(d)$ follows from them. It appears that one has to prove the following nontrivial

Identity: The coefficient of $x^n$ in $Td(x)^{n+1}e^{dx}$ equals
$$frac{1}{n!}
Td(partial /partial h_0) dots Td(partial /partial h_n)
(d+h_0+dots + h_n)^n |_{h_0=dots h_n=0}$$

A complete answer to my question would include a proof of this identity or a reference to where this is shown. (I did not find it in the cited papers.) I removed the acceptance to encourage a more complete explanation.