The Weierstrass Approximation Theorem states that any continuous function can be approximated on a bounded interval $[a, b]$ as closely as desired by a polynomial. In practice, instead of using Taylor series expansions, one uses polynomial minimax approximations $p(x)$ which approximate the function on the chosen interval with the minimal maximum error (thus the name).
When one looks at the graph of the error function $f(x) - p(x)$ one observes that it oscillates such that all extrema are of equal magnitude and neighboring extrema have opposite sign (equi-oscillation property), and that there are a total of $n+2$ extrema for a polynomial of degree $n$. This behavior of the error function can be used to construct the minimax polynomial, as also noted in gnasher729's answer.
The standard algorithm used to do this is called the Remez algorithm, named after the Russian mathematician who published it in 1927. This functionality can be found in ready-to-use form in widely used tools such as Maple and Mathematica, as well as the software tool Sollya, which provides an fpminimax command to generate minimax polynomial approximations. One can also create one's own implementation of the Remez algorithm, and I have used such code for this answer.
While the minimax polynomial is well defined mathematically, when the coefficients need to be stored in a finite-precision floating-point format and the polynomial evaluated is evaluated in finite-precision floating-point arithmetic, the error
function $(f(x) - p(x)$ becomes deformed and no longer satisfies the equi-oscillation criterion. That means we wind up with a polynomial that is merely close to the true minimax polynomial. Many tools do not adjust for this effect at all, delivering coefficients based on high-precision internal computations. Sollya takes the finite-precision storage of the coefficients into consideration, and thus can often deliver superior results.
A minimax polynomial of degree five suffices to approximate $x^{2.2}$ to within the accuracy desired by OP. $f(x)=p(x)+\epsilon, \:\: |\epsilon| < 7.4\cdot10^{-5}$. We can observe the equi-oscillation property of the error function $f(x) - p(x)$ and its seven extrema by plotting it, e.g. with Wolfram Alpha. Exemplary C code looks as follows:
/* Approximate x**(2.2) on [0,1] with minimax polynomial. Max. err < 7.4e-5 */
float pow_2p2_minimax (float x)
{
const float c0 = 1.06425546e-1f; // 0x1.b3eb46p-4
const float c1 = -3.56089801e-1f; // -0x1.6ca2cep-2
const float c2 = 5.86735249e-1f; // 0x1.2c6890p-1
const float c3 = 6.73461437e-1f; // 0x1.58cff0p-1
const float c4 = -1.05324369e-2f; // -0x1.59207cp-7
const float c5 = 7.38649687e-5f; // 0x1.35cfe8p-14
return ((((c0 * x + c1) * x + c2) * x + c3) * x + c4) * x + c5;
}
Note that the minimax polynomial does not return zero at the origin. If we desire this property, we can add this additional constraint by using $f(x) = x^{2}p(x) +\epsilon$. An example, with $|\epsilon| < 1.33 \cdot 10^{-4}$, is shown in this exemplary C implementation (plot of error function):
/* Approximate x**(2.2) on [0,1] with constrained minimax. Max. err < 1.33e-4 */
float pow_2p2_constraint (float x)
{
const float c0 = -3.66555989e-1f; // -0x1.775a74p-2
const float c1 = 1.19028902e+0f; // 0x1.30b6c8p+0
const float c2 = -1.55231142e+0f; // -0x1.8d6448p+0
const float c3 = 1.19035530e+0f; // 0x1.30bb20p+0
const float c4 = 5.38091123e-1f; // 0x1.1380aep-1
return x * x * ((((c0 * x + c1) * x + c2) * x + c3) * x + c4);
}
The computation of $x^{1/2.2}$ is a bit trickier. The easiest approach that does not suffer from premature underflow seems to be the following: Compute $y = \sqrt{\sqrt{x}} = x^{\frac{1}{4}}$, then compute $y^{\frac{20}{11}} = p(y) + \epsilon$ = $x^{\frac{5}{11}} = x^{1/2.2}$. The square root function is very fast on many modern processors, but slow on others, so this may not be a universal win, depending on how pow() is implemented. To achieve the error bound desired by OP, a minimax polynomial of degree six is needed, resulting in $|\epsilon| < 7.5 \cdot 10^{-5}$. Exemplary C code looks like so:
/* Approximate x**(1/2.2) on [0,1] with minimax polynomial. Max. err < 7.5e-5 */
float pow_0p4545_minimax (float x)
{
const float c0 = 3.07896197e-1f; // 0x1.3b4924p-2
const float c1 = -1.06079876e+0f; // -0x1.0f9082p+0
const float c2 = 1.48817670e+0f; // 0x1.7cf926p+0
const float c3 = -1.18180847e+0f; // -0x1.2e8b00p+0
const float c4 = 1.42678642e+0f; // 0x1.6d41e0p+0
const float c5 = 1.98969673e-2f; // 0x1.45fdeep-6
const float c6 = -7.45610159e-5f; // -0x1.38bb48p-14
/* compute x**0.25 */
x = sqrtf (sqrtf (x));
/* compute (x**0.25)**(20/11) = x**(1/2.2) */
return (((((c0 * x + c1) * x + c2) * x + c3) * x + c4) * x + c5) * x + c6;
}
