How to compute exp(x) for x in [0,ln(2)] in double precision with high accuracy, possibly at most 0.5 ulp error?

Step 1
I have found one nice formula for the exp function:
$$\begin{array}{}\text{(1)}& G(x)& =\frac{\coth <!--\; \; -->(\sqrt{x}/2)}{\sqrt{x}}-<!--\; -\; -->\frac{2}{x}\text{(2)}& c& =x-<!--\; -\; -->G(x^2)\times <!--\; \times \; -->x^2\text{(3)}& \exp <!--\; \; -->(x)& =1+(\frac{xc}{2-<!--\; -\; -->c}+x)\end{array}$$
This formula was used for interval [-ln(2)/2,ln(2)/2] but this is also good (transformation from this interval to [0,ln(2)] can be done without any errors. Unfortunately I don't have any reference on this scheme. I found this scheme using "reverse engeneering" of a library source.
The function G(x) is nearly linear in this interval. This function can be accurately approximated by a polynomial of order 4 using the Remez algorithm. Its Taylor expansion is
$$G(x)\approx <!--\; \approx \; -->\frac{1}{6}-<!--\; -\; -->\frac{x}{360}+\frac{x^2}{15120}-<!--\; -\; -->\frac{x^3}{604800}+\frac{x^4}{23950080}$$
Step 2
Condition number of this polynomial in given interval is very good, $1.012\times <!--\; \times \; -->{\mathrm{\gamma <!--\; \gamma \; -->}}_5$, ${\mathrm{\gamma <!--\; \gamma \; -->}}_N=2N/(1-<!--\; -\; -->2Nu)$.
Expression (2) cannot be computed accurately. However (3) suppresses most of numerical errors. For example condition number of (3) with respect to c is 0.0619, thus error if computing c, ${\mathrm{\delta <!--\; \delta \; -->}}_c$, contributes to the total error as $0.0619\times <!--\; \times \; -->{\mathrm{\delta <!--\; \delta \; -->}}_c$, similarly $xc/(2-<!--\; -\; -->c)$ is inaccurate, but error in computing this subexpression ${\mathrm{\delta <!--\; \delta \; -->}}_{xc}$ contributes to error of $xc/(2-<!--\; -\; -->c)+x$ approximately as $0.2\times <!--\; \times \; -->{\mathrm{\delta <!--\; \delta \; -->}}_{xc}$, which again is suppressed in the final sum (3) (approximately twice).
Final error of computing exp by (1)-(3) is less than 1 (maximum error found is 0.83).
It is possible to compute exp(x) even with 0.5 ulp error using the compensated Horner's scheme. See "Faithful Polynomial Evaluation with Compensated Horner Algorithm", P. Langlois, N. Louvet, 2006, for details. This scheme must however be corrected in order to take into account errors in computed coefficients of approximation polynomial. As the polynomial approximation of exp(x) we can use the Taylor expansion of order 27. Probably we can use lower order approximation, in tests I was not able to find a case with higher than 0.5 ulp even in case of approximation of order 20. Moreover this scheme is quite fast, at most 8-times slower, than library implementation of exp.