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>