SSE 2: Extremely slow floating-point `trunc`/`floor`/`ceil`/`round_ties_even`/`fract`

Nobody has claimed this yet.

Currently, SSE2 defers to the standard library (or libm) to implement certain floating-point operations: `trunc`/`floor`/`ceil`/`round_ties_even`/`fract`.

https://github.com/linebender/fearless_simd/blob/bd9f7c7be6b4a2c009558fffd5db67c54ff744f5/fearless_simd/src/generated/sse2.rs#L604-L647

While this produces correct results, it is also ***extremely*** slow. While SSE4.2 and up operate on 4 floats with a single instruction, SSE2 will perform one function call per f32. This can easily be ***100x slower*** for code using `f32x4` and wider.

A performance drop is to be expected when certain CPU features aren't available, but this cliff can absolutely kill performance.

## Micro benchmark

To show just how dramatic the difference is, here's a micro benchmark that just applies the operations onto 4096 floats using `f32x4::{fn}`. Floats are initialized to random bit patterns.

Benchmark code

```rs
use criterion::{criterion_group, criterion_main, Criterion};
use fearless_simd::{dispatch, f32x4, prelude::*, Level, Simd};
use rand::prelude::*;
use std::hint::black_box;

pub fn floating_ops_sse(c: &mut Criterion) {
let mut data = vec![0.0_f32; 4096];
for f in &mut data {
*f = f32::from_bits(rand::random::());
}

let level = Level::new();
let sse2 = level.as_sse2().unwrap().level();
let sse4_2 = level.as_sse4_2().unwrap().level();

macro_rules! bench_of {
($level:expr, $f:ident) => {
let name = format!("{} {}", stringify!($level), stringify!($f));
c.bench_function(&name, |b| {
let mut output = data.clone();

#[inline(always)]
fn process(simd: S, input: &[[f32; 4]], output: &mut [[f32; 4]]) {
for (i, o) in input.iter().zip(output.iter_mut()) {
*o = *f32x4::::simd_from(simd, *i).$f();
}
}

b.iter(|| {
let input = black_box(data.as_slice()).as_chunks::<4>().0;
let output = black_box(output.as_mut_slice()).as_chunks_mut::<4>().0;
dispatch!($level, simd => process(simd, input, output));
});
});
};
}

bench_of!(sse2, trunc);
bench_of!(sse4_2, trunc);
bench_of!(sse2, fract);
bench_of!(sse4_2, fract);
bench_of!(sse2, floor);
bench_of!(sse4_2, floor);
bench_of!(sse2, ceil);
bench_of!(sse4_2, ceil);
bench_of!(sse2, round_ties_even);
bench_of!(sse4_2, round_ties_even);
}

criterion_group!(benches, floating_ops_sse);
criterion_main!(benches);
```

Results on my Intel(R) Core(TM) i7-8700K CPU @ 3.70GHz CPU on Windows 10:

| Function | SSE 4.2 | SSE 2 | Perf drop |
| --- | ---: | ---: | :--- |
| `trunc` | 347.90 ns | 32.817 µs | 94x slower |
| `fract` | 387.02 ns | 33.925 µs | 88x slower |
| `floor` | 319.70 ns | 18.601 µs | 58x slower |
| `ceil` | 341.89 ns | 21.222 µs | 62x slower |
| `round_ties_even` | 319.66 ns | 133.41 µs | 417x slower |

Notes:
- All numbers are repeatable and consistent on my system. (Yes, `round_ties_even` is that slow.)
- SSE 4.2 versions are highly memory bound. Just loads and stores account for roughly 280 ns.

## Real-world example

In my use case, the hot loop uses `f32x16::trunc` and arithmetic operations in a roughly 1:8 ratio. Using `fearless_simd`'s `f32x16::trunc` made the entire loop 17x slower compared to the following fast approximation:

```rs
/// A fast version of `trunc` that works correctly for all x: -2^31 < x < 2^31.
/// Results are bit-exact for input values in the range of i32, and unspecified for input values outside of it.
#[inline(always)]
fn fast_trunc(x: f32x16) -> f32x16 {
#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
if matches!(x.simd.level(), Level::Sse2(_)) {
return i32x16::truncate_from(x).to_float();
}
x.trunc()
}
```

Using the approximation, the performance on SSE2 and SSE4.2 is within ~10%. A much more acceptable performance drop.

Details

This is the function processing one iteration in my loop (simplified).

```rs
#[simd]
fn block_closest_error_sq(&self, _: S, b: f32x16) -> f32 {
// factor1, factor2, add1, and add2 are all f32s and the only fields of Self
let blend = fast_trunc(b.mul_add(self.factor1, self.add1).max(0.0).min(7.0));
let closest = blend.mul_add(self.factor2, self.add2);
let error = b - closest;
(error * error).reduce_sum()
}
```

Here's the assembly it compiles to:

```nasm
; SSE4.2
<::block_closest_error_sq::__FearlessDispatch>::call::entry::, ::block_closest_error_sq::{closure#0}, f32>:
sub rsp, 88
movaps xmmword ptr [rsp + 64], xmm10
movaps xmmword ptr [rsp + 48], xmm9
movaps xmmword ptr [rsp + 32], xmm8
movaps xmmword ptr [rsp + 16], xmm7
movaps xmmword ptr [rsp], xmm6
movaps xmm3, xmmword ptr [rcx]
movaps xmm2, xmmword ptr [rcx + 16]
movaps xmm1, xmmword ptr [rcx + 32]
movaps xmm0, xmmword ptr [rcx + 48]
movss xmm5, dword ptr [rdx]
movss xmm6, dword ptr [rdx + 4]
movss xmm4, dword ptr [rdx + 8]
movss xmm8, dword ptr [rdx + 12]
shufps xmm6, xmm6, 0
shufps xmm8, xmm8, 0
movaps xmm7, xmm3
mulps xmm7, xmm6
addps xmm7, xmm8
movaps xmm9, xmm2
mulps xmm9, xmm6
addps xmm9, xmm8
movaps xmm10, xmm1
mulps xmm10, xmm6
addps xmm10, xmm8
mulps xmm6, xmm0
addps xmm6, xmm8
xorps xmm8, xmm8
maxps xmm7, xmm8
maxps xmm9, xmm8
maxps xmm10, xmm8
maxps xmm6, xmm8
movaps xmm8, xmmword ptr [rip + __xmm@40e0000040e0000040e0000040e00000]
minps xmm7, xmm8
minps xmm9, xmm8
minps xmm10, xmm8
minps xmm6, xmm8
roundps xmm7, xmm7, 11
roundps xmm8, xmm9, 11
roundps xmm9, xmm10, 11
roundps xmm6, xmm6, 11
shufps xmm4, xmm4, 0
shufps xmm5, xmm5, 0
mulps xmm7, xmm4
addps xmm7, xmm5
mulps xmm8, xmm4
addps xmm8, xmm5
mulps xmm9, xmm4
addps xmm9, xmm5
mulps xmm4, xmm6
addps xmm4, xmm5
subps xmm3, xmm7
subps xmm2, xmm8
subps xmm1, xmm9
subps xmm0, xmm4
mulps xmm3, xmm3
mulps xmm2, xmm2
mulps xmm1, xmm1
addps xmm1, xmm3
mulps xmm0, xmm0
addps xmm0, xmm2
addps xmm0, xmm1
movshdup xmm1, xmm0
addps xmm1, xmm0
movaps xmm0, xmm1
unpckhpd xmm0, xmm1
addss xmm0, xmm1
movaps xmm6, xmmword ptr [rsp]
movaps xmm7, xmmword ptr [rsp + 16]
movaps xmm8, xmmword ptr [rsp + 32]
movaps xmm9, xmmword ptr [rsp + 48]
movaps xmm10, xmmword ptr [rsp + 64]
add rsp, 88
ret
```
```nasm
; SSE2 with fast_trunc approximation
::block_closest_error_sq:::
sub rsp, 72
movaps xmmword ptr [rsp + 48], xmm9
movaps xmmword ptr [rsp + 32], xmm8
movaps xmmword ptr [rsp + 16], xmm7
movaps xmmword ptr [rsp], xmm6
movups xmm4, xmmword ptr [rcx]
movaps xmm3, xmmword ptr [rdx]
movaps xmm2, xmmword ptr [rdx + 16]
movaps xmm1, xmmword ptr [rdx + 32]
movaps xmm0, xmmword ptr [rdx + 48]
movaps xmm5, xmm4
shufps xmm5, xmm4, 85
movaps xmm7, xmm4
shufps xmm7, xmm4, 255
movaps xmm6, xmm5
mulps xmm6, xmm3
addps xmm6, xmm7
movaps xmm8, xmm5
mulps xmm8, xmm2
addps xmm8, xmm7
movaps xmm9, xmm5
mulps xmm9, xmm1
addps xmm9, xmm7
mulps xmm5, xmm0
addps xmm5, xmm7
xorps xmm7, xmm7
maxps xmm6, xmm7
maxps xmm8, xmm7
maxps xmm9, xmm7
maxps xmm5, xmm7
movaps xmm7, xmmword ptr [rip + __xmm@40e0000040e0000040e0000040e00000]
minps xmm6, xmm7
minps xmm8, xmm7
minps xmm9, xmm7
minps xmm5, xmm7
cvttps2dq xmm6, xmm6
cvttps2dq xmm7, xmm8
cvttps2dq xmm8, xmm9
cvttps2dq xmm5, xmm5
cvtdq2ps xmm6, xmm6
cvtdq2ps xmm7, xmm7
cvtdq2ps xmm8, xmm8
cvtdq2ps xmm5, xmm5
movaps xmm9, xmm4
shufps xmm9, xmm4, 170
mulps xmm6, xmm9
addps xmm6, xmm4
mulps xmm7, xmm9
addps xmm7, xmm4
mulps xmm8, xmm9
addps xmm8, xmm4
mulps xmm9, xmm5
addps xmm9, xmm4
subps xmm3, xmm6
subps xmm2, xmm7
subps xmm1, xmm8
subps xmm0, xmm9
mulps xmm3, xmm3
mulps xmm2, xmm2
mulps xmm1, xmm1
addps xmm1, xmm3
mulps xmm0, xmm0
addps xmm0, xmm2
addps xmm0, xmm1
movshdup xmm1, xmm0
addps xmm1, xmm0
movaps xmm0, xmm1
unpckhpd xmm0, xmm1
addss xmm0, xmm1
movaps xmm6, xmmword ptr [rsp]
movaps xmm7, xmmword ptr [rsp + 16]
movaps xmm8, xmmword ptr [rsp + 32]
movaps xmm9, xmmword ptr [rsp + 48]
add rsp, 72
ret
```
```nasm
; SSE2 with f32x16::trunc
::block_closest_error_sq:::
sub rsp, 280
movaps xmmword ptr [rsp + 256], xmm15
movaps xmmword ptr [rsp + 240], xmm14
movaps xmmword ptr [rsp + 224], xmm13
movaps xmmword ptr [rsp + 208], xmm12
movaps xmmword ptr [rsp + 192], xmm11
movaps xmmword ptr [rsp + 176], xmm10
movaps xmmword ptr [rsp + 160], xmm9
movaps xmmword ptr [rsp + 144], xmm8
movaps xmmword ptr [rsp + 128], xmm7
movaps xmmword ptr [rsp + 112], xmm6
movups xmm1, xmmword ptr [rcx]
movaps xmmword ptr [rsp + 32], xmm1
movaps xmm5, xmmword ptr [rdx]
movaps xmmword ptr [rsp + 96], xmm5
movaps xmm4, xmmword ptr [rdx + 16]
movaps xmmword ptr [rsp + 80], xmm4
movaps xmm3, xmmword ptr [rdx + 32]
movaps xmmword ptr [rsp + 64], xmm3
movaps xmm2, xmmword ptr [rdx + 48]
movaps xmmword ptr [rsp + 48], xmm2
movaps xmm6, xmm1
shufps xmm6, xmm1, 85
shufps xmm1, xmm1, 255
movaps xmm0, xmm6
mulps xmm0, xmm5
addps xmm0, xmm1
movaps xmm9, xmm6
mulps xmm9, xmm4
addps xmm9, xmm1
movaps xmm8, xmm6
mulps xmm8, xmm3
addps xmm8, xmm1
mulps xmm6, xmm2
addps xmm6, xmm1
xorps xmm1, xmm1
maxps xmm0, xmm1
maxps xmm9, xmm1
maxps xmm8, xmm1
maxps xmm6, xmm1
movaps xmm1, xmmword ptr [rip + __xmm@40e0000040e0000040e0000040e00000]
minps xmm0, xmm1
minps xmm9, xmm1
minps xmm8, xmm1
minps xmm6, xmm1
movshdup xmm13, xmm0
movaps xmm14, xmm0
unpckhpd xmm14, xmm0
movaps xmm15, xmm0
shufps xmm15, xmm0, 255
movshdup xmm12, xmm9
movaps xmm11, xmm9
unpckhpd xmm11, xmm9
movaps xmm10, xmm9
shufps xmm10, xmm9, 255
call truncf
movaps xmm7, xmm0
movaps xmm0, xmm13
call truncf
movaps xmm13, xmm0
movaps xmm0, xmm14
call truncf
movaps xmm14, xmm0
movaps xmm0, xmm15
call truncf
unpcklps xmm14, xmm0
unpcklps xmm7, xmm13
movlhps xmm7, xmm14
movaps xmm0, xmm9
call truncf
movaps xmm9, xmm0
movaps xmm0, xmm12
call truncf
movaps xmm12, xmm0
movaps xmm0, xmm11
call truncf
movaps xmm11, xmm0
movaps xmm0, xmm10
call truncf
unpcklps xmm11, xmm0
unpcklps xmm9, xmm12
movlhps xmm9, xmm11
movshdup xmm13, xmm8
movaps xmm14, xmm8
unpckhpd xmm14, xmm8
movaps xmm15, xmm8
shufps xmm15, xmm8, 255
movshdup xmm12, xmm6
movaps xmm11, xmm6
unpckhpd xmm11, xmm6
movaps xmm10, xmm6
shufps xmm10, xmm6, 255
movaps xmm0, xmm8
call truncf
movaps xmm8, xmm0
movaps xmm0, xmm13
call truncf
movaps xmm13, xmm0
movaps xmm0, xmm14
call truncf
movaps xmm14, xmm0
movaps xmm0, xmm15
call truncf
unpcklps xmm14, xmm0
unpcklps xmm8, xmm13
movlhps xmm8, xmm14
movaps xmm0, xmm6
call truncf
movaps xmm6, xmm0
movaps xmm0, xmm12
call truncf
movaps xmm12, xmm0
movaps xmm0, xmm11
call truncf
movaps xmm11, xmm0
movaps xmm0, xmm10
call truncf
unpcklps xmm11, xmm0
unpcklps xmm6, xmm12
movlhps xmm6, xmm11
movaps xmm1, xmmword ptr [rsp + 32]
movaps xmm0, xmm1
shufps xmm0, xmm1, 170
shufps xmm1, xmm1, 0
mulps xmm7, xmm0
addps xmm7, xmm1
mulps xmm9, xmm0
addps xmm9, xmm1
mulps xmm8, xmm0
addps xmm8, xmm1
mulps xmm0, xmm6
addps xmm0, xmm1
movaps xmm2, xmmword ptr [rsp + 96]
subps xmm2, xmm7
movaps xmm1, xmmword ptr [rsp + 80]
subps xmm1, xmm9
movaps xmm5, xmmword ptr [rsp + 64]
subps xmm5, xmm8
movaps xmm4, xmmword ptr [rsp + 48]
subps xmm4, xmm0
mulps xmm2, xmm2
movaps xmm3, xmm2
mulps xmm1, xmm1
movaps xmm2, xmm1
movaps xmm0, xmm5
mulps xmm0, xmm5
addps xmm0, xmm3
movaps xmm1, xmm0
movaps xmm0, xmm4
mulps xmm0, xmm4
addps xmm0, xmm2
addps xmm0, xmm1
movshdup xmm1, xmm0
addps xmm1, xmm0
movaps xmm0, xmm1
unpckhpd xmm0, xmm1
addss xmm0, xmm1
movaps xmm6, xmmword ptr [rsp + 112]
movaps xmm7, xmmword ptr [rsp + 128]
movaps xmm8, xmmword ptr [rsp + 144]
movaps xmm9, xmmword ptr [rsp + 160]
movaps xmm10, xmmword ptr [rsp + 176]
movaps xmm11, xmmword ptr [rsp + 192]
movaps xmm12, xmmword ptr [rsp + 208]
movaps xmm13, xmmword ptr [rsp + 224]
movaps xmm14, xmmword ptr [rsp + 240]
movaps xmm15, xmmword ptr [rsp + 256]
add rsp, 280
ret
```

## Suggestion

Either provide faster approximations for these operations or implement the exact methods more efficiently.

For `trunc`, I believe my approximation could be generalized to all f32s. The basic idea is that all f32s with an absolute value >=2^23 are integers, so we can return the f32 as is if the exponent is >=150 (which are all non-integer values + infinities + NaN). However, I'm unsure what the semantics of `f32::trunc` are for quiet and signaling NaNs, so that might complicate things.

For `round_ties_even`, I believe the old trick `(x_f32 + 8388608.0) - 8388608.0` could emulate it. Sign and values `x > 8388608` have to be handled carefully. We also have to guard against `fast-math` and x86 FPCR with different rounding modes.

I believe `floor` and `ceil` could be implemented in terms of `trunc`.

Start with the SSE2 implementations in fearless_simd/src/generated/sse2.rs around lines 604-647, then compare them with the SSE4.2 implementations and the supplied fast_trunc approximation. Use the provided Criterion benchmark over 4096 f32 values to verify that trunc, floor, ceil, round_ties_even, and fract remain correct while avoiding the severe SSE2 performance cliff.

Written by the indexing model from the issue text.