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