It would be nice if we knew the error bounds for each of the
approximation methods. Do we know how the coefficients were generated?
When switching from one method to another, it might be nice to have a
range where the results slowly blend from one method to another:
if x < a: return f1(x)
if x > b: return f2(x)
t = (x - a) / (b - a)
return (1-t)*f1(x) + t*f2(x)
This minimizes discontinuities in the first and second derivatives.
The lowword/highword macros look to be tightly tied to the internal
processor representation for floats. It would be more portable and
maintainable to replace that bit twiddling with something based on frexp
().