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

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

• If you are looking for an implementation of the Langevin function, check out my question here Jan 19 '20 at 6:11
• @njuffa. I did not know that this function was already well-studied; for me it came up in the derivative of $\log( (\exp(x)-1) / x )$. Thank you for identifying it! In your exemplary C function (with a different minimax polynomial near 0) the switch-over is 1.8125. Can you give some insight on a general procedure for choosing the optimal switch-over ? Jan 19 '20 at 17:12
• I determined the switchover point experimentally, trying to achieve similar amounts of ulp error on either side. I do not recall the specifics of this case, but I often use an iterative process as part of which the switchover point shifts as I create increasingly accurate approximations for two or more approximation intervals. I also try to make switchover points small fractions with power-of-two denominators, but that is more of a personal quirk. Jan 19 '20 at 19:45
• @njuffa. The only way I can think of to determine the ulp error is to compare the computed values with some type of "gold standard". Where does this "gold standard" for the Langevin function come from ? If there is another way to determine the ulp error, what is it ? Jan 21 '20 at 4:25
• I used an arbitrary-precision library for my golden test reference. This allows the use of naive mathematical formulas in the presence of subtractive cancellation. For arguments $\lt 10^{-10}$ I used the approximation $x/3$. Jan 21 '20 at 4:41

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