As far as I understand, there is no point in considering approximations that are slower than the standard acos() or acosf() functions. Achieving same performance for correctly rounded double-precision values ​​is extremely difficult, if not impossible, but it is quite possible to improve performance for values with error ​​close to one of single-precision format. Therefore, even those approximations that seem successful should be tested for performance.

Since the arccosine of x has an unbounded derivative at points x=+-1, the approximated function should be transformed so that it becomes sufficiently smooth. I propose to do this as follows (I think this is not a new way): is constructed an approximation of the function

f(t) = arccos(t^2)/(1-t^2)^0.5

using the Padé-Chebyshev method, where t=|x|^0.5, -1<=t<=1. The function f(t) is even, fairly smooth, and can be well approximated by both polynomial and a fractional rational functions. The approximation is as follows:

f(t) ≈ (a0+a1*t^2+a2*t^4+a3*t^6)/(b0+b1*t^2+b2*t^4+b3*t^6) = p(t)/q(t).

Considering the relationship between the variables t and x, we can write:

f(x) ≈ (a0+a1*|x|+a2*|x|^2+a3*|x|^3)/(b0+b1*|x|+b2*|x|^2+b3*|x|^3) = p(x)/q(x).

After calculating the function f(x), the final result is obtained using one of the formulas:

arccos(x) = f(x)*(1-|x|)^0.5 at x>=0;

arccos(x) = pi-f(x)*(1-|x|)^0.5 at x<=0.

The coefficients of the fractional rational function f(x), providing a maximum relative error of 8.6E-10, are follows:

a0 = 1.171233654022217, a1 = 1.301361441612244, a2 = 0.3297972381114960, a3 = 0.01141332555562258;

b0 = 0.7456305027008057, b1 = 0.9303402304649353, b2 = 0.2947896122932434, b3 = 0.01890071667730808.

These coefficients are specially selected for calculations in single precision format.

An example of code implementation using the proposed method can be found in the adjacent topic Fast Arc Cos algorithm?