Euler-Maclaurin Sum Formula
---------------------------
The Euler-Maclaurin Sum Formula is based on the fact that, formally, the
Taylor series for F(x-1) is

     -D
    e   F(x)                                                                [1]

where D is d/dx.  Thus, if f(x) = F(x) - F(x-1), then

                 -D
    f(x) = (1 - e  ) F(x)                                                   [2]

If we "solve" [2] for F(x), we get

               D     -1
    F(x) =  ------- D  f(x)                                                 [3]
                 -D
            1 - e

where D^-1 is indefinite integration (complete with constant of integration).
The differential operator D/(1 - e^-D) is

        1      1  2    1   4     1    6      1     8       1     10
    1 + - D + -- D  - --- D  + ----- D  - ------- D  + -------- D   - ...   [4]
        2     12      720      30240      1209600      47900160

We get [4] by plugging D into x/(1 - e^-x).  Thus, asymptotically,

     n
    ---           -1       1         1          1   3         1    5
    >   f(k) ~   D  f(n) + - f(n) + -- Df(n) - --- D f(n) + ----- D f(n) - ...
    ---                    2        12         720          30240
    k=1

To compute more terms in [4],

                                     oo
       x      x   x       x     x   ---     2n
    ------- = - + - coth( - ) = - + >   a  x
         -x   2   2       2     2   ---  n
    1 - e                           n=0

                  1
    a  = 1, a  = --
     0       1   12

              n-1
          -1  ---
    a  = ---- >   a  a
     n   2n+1 ---  k  n-k
              k=1

Shifted Euler-Maclaurin Sum Formula
-----------------------------------
Start with equation [3] and "multiply" by e^{D/2}/e^{D/2}:

                  D       -1  D/2
    F(x) =  ------------ D   e    f(x)                                      [5]
             D/2    -D/2
            e    - e

where e^{D/2} f(x) = f(x+1/2) and D^-1 is indefinite integration.
The differential operator D/(e^{D/2}-e^{-D/2}) is

         1  2    7    4     31    6      127     8       73      10
    1 - -- D  + ---- D  - ------ D  + --------- D  - ---------- D   + ...   [6]
        24      5760      967680      154828800      3503554560

We get [6] by plugging D into x/2 csch(x/2).  Thus, asymptotically,

     n
    ---           -1            1              7    3
    >   f(k) ~   D  f(n+1/2) - -- Df(n+1/2) + ---- D f(n+1/2) - ...
    ---                        24             5760
    k=1

To compute more terms in [6],

                  oo
    x       x     ---     2n
    - csch( - ) = >   a  x
    2       2     ---  n
                  n=0

    a  = 1
     0

            n
           ---           1
    a  = - >   a    -----------
     n     ---  n-k 4^k (2k+1)!
           k=1

Euler-Maclaurin Sum Formula with Remainder
------------------------------------------
If f(x) = F(x) - F(x-1), then the Taylor series with remainder says

            n                                   n
           ---     k-1 1   k        |\x  (x-1-t)   n+1
    f(x) = >   (-1)    -- D F(x) +  |    -------  D   F(t) dt               
           ---         k!          \|x-1    n!
           k=1

Taking m derivatives of both sides, m <= n, yields

             n-m                                    n-m
     m       ---     k-1 1   k+m        |\x  (x-1-t)     n+1
    D f(x) = >   (-1)    -- D   F(x) +  |    -------    D   F(t) dt
             ---         k!            \|x-1  (n-m)!
             k=1

One term from the summation gets cancelled by the integral with each
derivative.

Notice that the differential operator inside the sum for f(x) is the
first n terms of 1 - e^{-D} (degrees 1 to n).  Therefore, if we apply
E_n(D), the first n terms of E(D) = D/(1 - e^{-D}) (degrees 0 to n-1),
to f(x), we get
                                              n
                          |\x          (x-1-t)    n+1
    E (D) f(x) = DF(x) +  |    ( E (D) ------- ) D   F(t) dt
     n                   \|x-1    n       n!

With a change of variables, an antiderivative, and a reverse change of
variables, we have
                                                n
           -1               |\x          (x-1-t)    n
    E (D) D  f(x) = F(x) +  |    ( E (D) ------- ) D F(t) dt
     n                     \|x-1    n       n!

Let B_n(x) = E_n(D) x^n/n!, then the formula with error term becomes

           -1               |\x             n
    E (D) D  f(x) = F(x) +  |    B (x-1-t) D F(t) dt
     n                     \|x-1  n

Some properties of E(D), E_n(D), and B_n(x) are derived in the next
section.

Methods of estimating the size of the error term can vary depending on
the functions involved.  Usually, high derivatives of f are decreasing
functions, so we can use the fact that f(x) = F(x) - F(x-1) = DF(y) for
some y between x-1 and x to say that for t between x-1 and x,

       n             n-1
    | D F(t) | <= | D   f(x-1) | 

Let

    M  =  max  |B (t)|
     n   [-1,0]  n

Then the error term is bounded by

          n-1
    M  | D   f(x-1) |
     n

The greatest order of D applied to f in E_n(D) D^{-1} f(x) is n-2.
Usually, D^{n-1} f(x) decreases faster than D^{n-2} f(x), so the error
term decreases faster than the last term used in the Euler-Maclaurin
Sum Formula.

Properties of E(D), E_n(D), and B_n(x)
--------------------------------------
Notice that

           D   e^{D/2} + e^{-D/2} D      D/2
    E(D) - - = ------------------ - = ---------
           2   e^{D/2} - e^{-D/2} 2   tanh(D/2)

and

          -D/2           D               D/2
    E(D) e     = ------------------ = ---------
                 e^{D/2} - e^{-D/2}   sinh(D/2)

Therefore, both E(D) - D/2 and E(D) e^{-D/2} are even operators.  We
will not worry about convergence with infinite series in D since they
will only operate on polynomials, where convergence is not an issue.
In all other cases, we will work with truncated series (polynomials)
in D.

By definition, E(D) - E_n(D) consists only of terms of order at least
n.  Furthermore, since E(D) - D/2 is even, for n > 1, E(D) - E_n(D) is
an even operator.

Because E(D) - E_n(D) consists only of terms of order at least n,
(E(D) - E_n(D)) x^n/n! must be constant.  Furthermore, if n is odd and
not 1, (E(D) - E_n(D)) x^n/n! must be an odd function, and the only
odd, constant function is 0.  Therefore,

                 x^n                    x^n
    B (x) = E(D) --- - ( E(D) - E (D) ) ---
     n            n!             n       n!

                 x^n
          = E(D) --- - c
                  n!    n

where c_n = 0 when n is odd and not 1.  Apply e^{-D/2} to both sides:

            1           -D/2 x^n
    B ( x - - ) = E(D) e     --- - c
     n      2                 n!    n

This says that B_n(x-1/2) is even if n is even and odd if n is odd and
not 1.

Since the greatest order of derivative in E_n(D) is n-1, we have that
each term in B_n(x) = E_n(D) x^n/n! has a factor of x, and therefore,
B_n(0) = 0.  When n is not 1, B_n(x-1/2) is either even or odd; thus,
we also have that B_n(-1) = 0, and if n is odd, B_n(-1/2) = 0.

As noted above, B_n(x) = E(D) x^n/n! - c_n.  Therefore,

                             -D         x^n
    B (x) - B (x-1) = ( 1 - e  ) ( E(D) --- - c )
     n       n                           n!    n

                          x^n
                      = D --- - 0
                           n!

                        x^{n-1}
                      = -------
                         (n-1}!

Since B_n(0) = 0, we get that for integer k,

                    k
              1    ---  n-1
    B (k) = ------ >   j
     n      (n-1)! ---
                   j=1