For approximation of some trigonometric functions you can evaluate an approximating polynomial which can be faster if you can live with restrictions in the parameter range. This article is not about acquiring the coefficients of these polynomials but rather the implementation of their evaluation.
For a fast sine approximation I use the following polynomial as an example (from 3D Game Engine Architecture, David H. Eberly, second printing, ISBN 978-0-12-229064-0):
sinx/x = 1 - 0.16605x^2 + 0.00761x^4 + e(x)
The input range is [0,pi/2] and the error: |e(x)| <= 1.7e-4
As there are a few of these approximations and typing out all of those multiplies manually is error-prone, I use recursive template instantiation to accomplish this task. This does not compromise runtime performance as you will see later.
The implementation of this particular polynomial looks like this:
static inline __m128 fast_sin_0(__m128 angles) {
static const float coeffs[] = { -0.16605f, 0.00761f };
return angles * (
ones() +
polynomialeval_even_packed_t<
sizeof(coeffs)/sizeof(coeffs[0]),
float,__m128
>()(angles, coeffs)
);
}
The polynomialeval_even_packed_t is a templated struct implementing the call operator:
template<unsigned N, typename real_t, typename packed_t>
struct polynomialeval_even_packed_t {
inline packed_t operator()(packed_t x, const real_t * coeffs) {
return coeffs[N-1] * expeval_packed_t<N*2,real_t,packed_t>()(x) +
polynomialeval_even_packed_t<N-1,real_t,packed_t>()(x, coeffs);
}
};
template<typename real_t, typename packed_t>
struct polynomialeval_even_packed_t<0,real_t,packed_t> {
inline packed_t operator()(packed_t x, const real_t * coeffs) {
return math_t<real_t,packed_t>::zeroes();
}
};
This struct returns sum[0<=i<N](coeffs[i]*x^(n*2+2)), which means for the case of fast_sin_0() which uses a coefficient array of size 2:
sum[0<=i<2](coeffs[i]*x^(n*2+2)) = coeffs[0]*x^2 + coeffs[1]*x^4
The expeval_packed_t struct returns x^N.
Now lets look at the code generated by VC++12 for x64:
?fast_sin_0@?... PROC 00000 0f 28 21 movaps xmm4, XMMWORD PTR [rcx] 00003 0f 57 c9 xorps xmm1, xmm1 00006 0f 28 c4 movaps xmm0, xmm4 00009 0f 59 c4 mulps xmm0, xmm4 0000c 0f 28 d8 movaps xmm3, xmm0 0000f 0f 59 05 00 00 00 00 mulps xmm0, XMMWORD PTR __xmm@be2a0903be2a0903be2a0903be2a0903 00016 0f 59 dc mulps xmm3, xmm4 00019 0f 58 c1 addps xmm0, xmm1 0001c 0f 59 dc mulps xmm3, xmm4 0001f 0f 59 1d 00 00 00 00 mulps xmm3, XMMWORD PTR __xmm@3bf95d4f3bf95d4f3bf95d4f3bf95d4f 00026 0f 58 c3 addps xmm0, xmm3 00029 0f 58 05 00 00 00 00 addps xmm0, XMMWORD PTR ?ones_@?1??... 00030 0f 59 c4 mulps xmm0, xmm4 00033 c3 ret 0 ?fast_sin_0@?... ENDP
This is not bad and it was -- as requested -- inlined (not shown here). The conversion of the scalar coefficients to packed values was done at compile time. There's an unneeded addition of zero involved which is used to terminate the template instantiation recursion. I compiled with /fp:precise, so compiling with /fp:fast might lead to the removal of that addition.
