I'll share it. I think people will appreciate the simplicity and hopefully make use of the idea.
First natural thing to do: keep repeating integration by parts on the integral of 1*f.
Ignoring the error and letting M = infinity gives you infinitely many constants of integration (one for each p_m, m>0) to freely choose.
What would be a "convenient" choice for all of them?
Wouldn't it be awesome if they were chosen so that if the formula was applied to [0,1] and added to the formula applied to [1,2] then there would be maximal cancellation for all derivatives of f at 1? (Telescoping)
The euler maclaurin summation formula and the generating function for bernoulli polynomials/numbers are obvious from here.
I'll let you derive the "midpoint" method which has a much better error term.
The trick for the midpoint method is to break [0,1] in half and use two alternating methods for the interval [0,1/2].
P1 = k1(t)e^(xt) and P2 = k2(t)e^(xt).
P1(1/2,t)-P2(0,t) = t and P2(1/2,t) - P1(0,t) = 0.