Based on the double precision's sin/cos fast path algorithm:
Step 1: Perform range reduction `y = x mod pi/8` with target errors <
2^-54.
This is because the worst case mod pi/8 for single precision is ~2^-31,
so to have up to 1 ULP errors from
the range reduction, the targeted errors should `be 2^(-31 - 23) =
2^-54`.
Step 2: Polynomial approximation
We use degree-5 and degree-4 polynomials to approximate sin and cos of
the reduced angle respectively.
Step 3: Combine the results using trig identities
```math
\begin{align*}
\sin(x) &= \sin(y) \cdot \cos(k \cdot \frac{\pi}{8}) + \cos(y) \cdot \sin(k \cdot \frac{\pi}{8}) \\
\cos(x) &= \cos(y) \cdot \cos(k \cdot \frac{\pi}{8}) - \sin(y) \cdot \sin(k \cdot \frac{\pi}{8})
\end{align*}
```
Overall errors: <= 3 ULPs for default rounding modes (tested
exhaustively).
Current limitation: large range reduction requires FMA instructions for
binary32. This restriction will be removed in the followup PR.
---------
Co-authored-by: Petr Hosek <phosek@google.com>