evaluating $\coth(x) - 1/x$ for real $x$, on 2 "pieces"

Question

The function $\coth(x) - 1/x$ has a removable singularity at 0. Its Taylor series is: $$ \coth(x) - 1/x = \frac{x}{3} - \frac{x^3}{45} + \frac{2x^5}{945} + \ldots $$ I would like to evaluate the first 3 terms of the Taylor series for $|x| \le \delta$, and just do a straight calculation of $\coth(x)$ and $x^{-1}$ for $|x|>\delta$.

What is the numerically correct way to select $\delta$, so the error over some fixed finite interval, say [-10,10], is minimized ?

Assume $x$ is 64-bit IEEE floating point, and that $\coth(x)$ is taken from a high-quality 64-bit precision library.

If $\delta$ is too small, then there is cancellation error for $x$ just a bit larger than $\delta$. And if $\delta$ is too large, then there is approximation error for $x$ just a bit less than $\delta$.

Surely there is a smart way to analyze this without using trial-and-error.

Added later: Thanks to @njuffa, I have learned that this is the Langevin function and that a different polynomial, or perhaps a continued fraction expansion, near 0 might be more accurate. So a more general question is: after an approximate formula near 0 is chosen, what is a good procedure for finding the best switch-over point $\delta$ ?

Answer 1

Following @njuffa, I used the arbitrary precision R package Rmpfr to compute a golden test reference, with 1024 bits of precision. Here is a plot of relative error, for both the straightforward definition, and the Taylor series (quintic polynomial) approximation (both using double-precision):

The optimal $\delta$ is where the two curves cross, and is about $\delta \approx 0.03$. The relative error at $\delta$ is about $6.6 \times 10^{-13}$, which is about 3000 times the machine epsilon = 2.220446e-16. The machine epsilon appears in the plot as the dotted horizontal line. The plot shows that subtractive cancellation occurs for all $x$ less than $\approx 2$, which confirms what @njuffa wrote in his previous question. Unity is too small.

These 2 "curves" are simpler than I expected to be, and so it may be possible to derive simple upper bound functions for them from first principles. If so then one could solve for $\delta$ from first principles, and with no experimentation. This is what I had in mind when I asked the question. However, using the arbitrary-precision library is much easier and was actually fun.

I also tried out the continued fraction expansion of the Langevin function on Wikipedia given here. The denominators are consecutive odd integers starting at 3. For results comparable to @njuffa 's exemplary C function, one must go all the way out in the continued fraction to denominator 21 (not 25 I made a programming error), and use the switch-over point $\delta \approx 1.8$.