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^{2^{-100}} \approx 1 + 2^{-100}$ and then subtracting 1 results in a huge (relative) error.

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. (The IEEE floating point standard "recommends" to include it in its implementations.)

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