Our approach to computing numerical values of these sums involves the
compound application of the Euler-Maclaurin summation formula (see
[2, p. 289,]), which can be stated as follows. Suppose
has at least 2p + 2 continuous derivatives on . Let D be
the differentiation operator, let denote the k-th Bernoulli
number, and let denote the k-th Bernoulli polynomial.
Then
where the remainder is given [2, p. 289,] by