I suspect that using the expm1
function would give you a better result.
Computing $e^x -1$ with the trivial formula in machine precision gives you only limited accuracy for small inputs: the fundamental reason is that there are only "few" floating point numbers around 1, and when you first compute $e^x$ the machine has to approximate your result to the closest floating point number. For instance, the next number after $1$ is $1+2^{-52}$, so every number in $[1, 1+2^{-52}]$ has to be replaced with one of the two extremes. So computing, for instance $e^x-1$ for $x=2^{-100}$ with good accuracy is a lost cause.
Hence there is the need for a separate library function that computes $e^x-1$ natively and is accurate also for small inputs (it can be implemented using a Taylor expansion for small inputs, for instance). Hence expm1
, which is in most programming languages and is an optional part of the IEEE floating point format.
(Similarly, there is log1p
to compute $\log(1+x)$. These functions are very handy when working with very small probabilities in log space, for instance.)
In any case, if you know that your function has a multiple zero, there is no reason not to use a Newton variant designed for multiple zeros.