Is expm1 the right primitive?

Question

I'm writing some code to calculate $\int_0^1 e^{ax} \mathrm{d} x$. Annoyingly there does not seem to be a way of doing this without if statements:

def my_integral(a):
    if a == 0:
        return 1.
    else:
        return expm1(a)/a

(I know comparing to zero looks worrying, but this correctly returns 1 for values like 1e-50)

It got me thinking; expm1 is really standard, but why isn't the right primitive expm1_over(x) = expm1(x)/x?

expm1_over can do everything that expm1 can, and more.

Am I missing some drawbacks?

Answer 1

Yes, it is very well arguable that $\frac{e^x-1}{x}$ is ``as right as'' expm1. It appears frequently in applications; see for instance Sections 2.1 and 10.7.4 of Higham's Function of Matrices, which mention exponential integrators and matrix functions; this function is called $\psi_1(x)$ there. Your example, integrating exponentials, is another one, and it is very relevant in probability applications as well. Personally, I have encountered $\psi_1$ more often than expm1 in my research.

As @WolfangBangerth notes in his answer, it is not really a big deal in practice, because one can be computed from the other with high relative accuracy. It seems to me that the only annoyance is the one you point out.

Answer 2

expm1 is a primitive because for small $x$, the Taylor series for $e^x=1+x+x^2+\cdots$ adds small terms to one -- which implies that if $x$ is small enough, you lose almost all digits of accuracy other than the leading one and so the expression $e^x-1$ can not be accurately computed by first evaluating the exponential.

The fact that you happen to have an application that requires some other primitive has no impact on the usefulness of expm1. I will, however, point out that there is no reason to invent your expm1_over primitive because the operation (e^x-1)/x can be accurately computed using first the expm1 primitive followed by the division, without loss of accuracy.