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.
        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?

  • $\begingroup$ What's wrong with this question? $\endgroup$ – user357269 Nov 26 '19 at 16:04
  • $\begingroup$ This question has nothing to do with computational sciences. It's more about the implementation of a certain library, which of course is off-topic here. $\endgroup$ – Alone Programmer Nov 26 '19 at 16:09
  • 3
    $\begingroup$ Numerical stability is a pretty important component of computational science. That said, I'm also unclear what the question is. $\endgroup$ – Richard Nov 26 '19 at 16:23
  • 2
    $\begingroup$ In my opinion, there is nothing wrong with the question (I edited the tags though), it is certainly on-topic here, even with or without mentioning of numpy (which has no real importance to the question). $\endgroup$ – Anton Menshov Nov 26 '19 at 18:13
  • 1
    $\begingroup$ @AloneProgrammer the list of on-topic is here. This particular question would correspond to multiple points: programming that requires computational science experience, is of numerical nature and is high-level wrt numpy or an abstract numerical library implementation. $\endgroup$ – Anton Menshov Nov 26 '19 at 18:35

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.

| cite | improve this answer | |

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.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.