  Bernstein's theorem (also called Hausdorff--Bernstein--Widder theorem)
enables the integral representation of a completely monotonic function. We
introduce a finite completely monotonic function, which is a completely
monotonic function with a finite positive integral interval of the integral
representation. We consider the exponential sum approximation of a finite
completely monotonic function based on the Gaussian quadrature with a variable
transformation. If the variable transformation is analytic on an open Bernstein
ellipse, the maximum absolute error decreases at least geometrically with
respect to the number of exponential functions. The maximization of the
decreasing rate of the error bound can be achieved by using a variable
transformation represented by Jacobi's delta amplitude function (also called dn
function). The error curve is expanded by introducing basis functions, which
are eigenfunctions of a fourth order differential operator, satisfy
orthogonality conditions, and have the interlacing property of zeros by
Kellogg's theorem.
