Differences between SIMD and non-SIMD algorithm results

[qarmin]

Problems found by tests in PR - #160

Below are visible potential fixes for such problems(AI guessing)

# Backend-vs-scalar comparison test findings

A differential test suite (every SIMD backend compared against the crate's own
portable scalar fallback, `#[cfg(test)] mod backend_comparison_tests` in each
`src/{sse,avx2,avx512,neon,sve2,wasm32}/*.rs` file) found the issues below.
Every affected test is currently marked `#[ignore = "known bug: ..."]` so CI
stays green; remove the `#[ignore]` once a fix lands to confirm it.

## Alpha premultiply / unpremultiply (SSE + AVX2, u8)

- **`_mm_div_by_255_epi16`** (`src/sse/alpha_u8.rs:44`) and its AVX2
  counterpart in `src/avx2/alpha_u8.rs` (used by `Avx2PremultiplyExecutor*`)
  round divide-by-255 differently than the scalar reference. The scalar
  `div_by_255` (`src/alpha_handle_u8.rs:47-49`) computes
  `t = v + 128; (t + (t >> 8)) >> 8`. The SIMD version instead adds `127` on
  one term and right-shifts the *unbiased* `v` on the other - a different
  additive constant on each term. This produces off-by-one results for part
  of the `u8 * u8` product range (e.g. `v = 128` gives `0` instead of `1`).
  - Tests: `sse::alpha_u8::sse_premultiply_matches_scalar_reference`,
    `avx2::alpha_u8::avx2_premultiply_matches_scalar_reference`.
  - Fix: change the SIMD formula to add the rounding bias (`128`) before both
    the direct term and the shifted term, matching the scalar reference
    exactly, or replace with a lookup-table-free implementation proven
    equivalent to `(v + 128 + ((v + 128) >> 8)) >> 8`.

- **`sse_unpremultiply_row`** (`src/sse/alpha_u8.rs:54`) and
  `Avx2DisassociateAlpha` / `Avx2DisassociateAlphaFast`
  (`src/avx2/alpha_u8.rs:151`, `:240`) use the low-precision approximate
  reciprocal instructions `_mm_rcp_ps` / `_mm256_rcp_ps` instead of exact
  division. The scalar reference uses an exact (or table-based) division, so
  results differ by a few units for some alpha values.
  - Tests: `sse::alpha_u8::sse_unpremultiply_matches_scalar_reference`,
    `avx2::alpha_u8::avx2_unpremultiply_quality_matches_scalar_reference`,
    `avx2::alpha_u8::avx2_unpremultiply_speed_matches_scalar_reference`.
  - Fix: either do a real division (`_mm_div_ps`/`_mm256_div_ps`), or run one
    Newton-Raphson refinement step after `_mm_rcp_ps`/`_mm256_rcp_ps` to
    bring the reciprocal to full precision before multiplying.

- **AVX2 u16 unpremultiply** (`src/avx2/alpha_u16.rs`, masks built around
  `is_zero_alpha_mask*` / `_mm256_blendv_epi8`) zeroes the RGB channels when
  alpha is `0`. The scalar u16 unpremultiply path leaves RGB untouched when
  alpha is `0` (only the *u8* scalar path zeroes on `alpha == 0`), so the two
  bit depths currently disagree with each other, and AVX2 disagrees with its
  own scalar fallback.
  - Test: `avx2::alpha_u16::avx2_unpremultiply_matches_scalar_reference`.
  - Fix: decide the intended contract for `alpha == 0` once (zero RGB or
    preserve it) and make the u8 scalar, u16 scalar, and every SIMD backend
    agree with it consistently.

- The same alpha==0 inconsistency class likely also affects
  **`src/neon/alpha_u16.rs`** (`neon::alpha_u16::neon_unpremultiply_matches_scalar_reference`)
  and the NEON fp16 path in **`src/neon/alpha_f16_full.rs`**
  (`neon::alpha_f16_full::neon_fp16_unpremultiply_matches_scalar_reference`) -
  same root cause, needs the same fix applied consistently across all
  backends/bit-depths once decided.

## SSE f32→f16 conversion drops the sign bit

- **`_mm_cvtps_ph_fallback`** (`src/sse/f16_utils.rs:144`), the software f32→f16
  conversion SSE always uses (the hardware `f16c`-accelerated path is not
  selected here), computes the exponent/mantissa bits and packs them into
  the result but never extracts or ORs in the IEEE-754 sign bit. Every
  negative `f32` therefore round-trips through this function as its
  **positive** counterpart. Since resampling kernels like Lanczos3 have
  negative lobes, this corrupts real f16 resize output whenever a negative
  kernel tap is involved - not just a test artifact.
  - Tests: `sse::rgba_f16::sse_row_matches_scalar_reference`,
    `sse::rgba_f16::sse_rows_4_matches_scalar_reference`,
    `sse::rgb_f16::sse_row_matches_scalar_reference`,
    `sse::rgb_f16::sse_rows_4_matches_scalar_reference`,
    `sse::vertical_f16::sse_f16_vertical_matches_scalar_reference`
    (one root cause, five failing tests).
  - Fix: extract the sign bit (`x_bits & 0x8000_0000`) from the input `f32`
    and OR it into the correct position of the packed `f16` result before
    returning.

## FMA-dependent (±1) rounding divergences, u16 (SSE/AVX2) - found by CI, not locally

These four were **not** caught during initial local development - the dev
machine used always has AVX2/FMA available, so the "SSE-only" and "AVX2
without FMA" code paths were being exercised, but the *compiler* was still
free to contract multiply+add into FMA instructions elsewhere in the
compilation unit, which coincidentally matched the scalar reference. CI's
`Testing x86 SSE isolated` and `Testing x86 AVX2 isolated (no FMA)` jobs
explicitly compile with `-C target-feature=...,-fma` and run under
`qemu-x86_64 -cpu Nehalem` / `-cpu Haswell-noTSX,avx2=on,fma=off` specifically
to catch this class of bug, and did:

- `sse::alpha_u16::sse_unpremultiply_matches_scalar_reference`
- `avx2::rgb_u16::avx2_row_matches_scalar_reference` (the single-row variant;
  `avx2_rows_4` was already known-bad, see below)
- `avx2::rgba_u16::avx2_default_row_matches_scalar_reference` /
  `avx2_default_rows_4_matches_scalar_reference` - note this is the
  **"default" (non-FMA-flagged) variant**, not the `_fma` one; it turns out
  the "default" AVX2 kernels are *not* actually bit-exact either once FMA
  contraction is truly disabled at the compiler level, contrary to what
  local-only testing suggested.

All four reproduce identically under
`RUSTFLAGS="-C target-feature=+sse4.1,-avx,-avx2,-fma" CARGO_TARGET_X86_64_UNKNOWN_LINUX_GNU_RUNNER="qemu-x86_64 -cpu Nehalem" cargo test --target x86_64-unknown-linux-gnu ...`
and the AVX2 equivalent - use these locally to verify a fix without waiting
on CI. Likely the same underlying rounding-direction root cause as the next
section; worth investigating together.

## Small (±1) fixed-point rounding divergences, u16 (SSE/AVX2/NEON)

A consistent pattern across many `u16` fixed-point SIMD paths: output
differs from the scalar reference by exactly `1` in a small number of output
elements, reproducible specifically with **Lanczos3** (a kernel with negative
taps) rather than with `Bilinear`/`Nearest`. This suggests a shared rounding
behavior difference when accumulating negative products (e.g. arithmetic
right-shift rounding toward negative infinity in the SIMD path vs
round-half-up in the scalar path, or an intermediate truncation difference).
Worth investigating as one shared root cause rather than N separate bugs.

Affected (all bit-exact fixed-point paths that should never differ from the
scalar reference by design):

- `avx2::rgb_u16::avx2_rows_4_matches_scalar_reference`
- `avx2::plane_u16::avx2_rows_4_matches_scalar_reference`
- `avx2::rgba_u16::avx2_fma_row_matches_scalar_reference` /
  `avx2_fma_rows_4_matches_scalar_reference` (only the **FMA** variant - the
  "default" non-FMA AVX2 variant is bit-exact, which supports the "rounding
  behavior differs under FMA" theory above)
- `avx2::vertical_u16::avx2_default_vertical_matches_scalar_reference` /
  `avx2_fma_vertical_matches_scalar_reference` (the >12-bit "high bit depth"
  vertical path, which internally accumulates in `f32` rather than fixed-point
  - see `src/avx2/vertical_u16.rs`)
- `sse::rgba_u16::sse_row_matches_scalar_reference` /
  `sse_rows_4_matches_scalar_reference`
- `sse::vertical_u16::sse_vertical_matches_scalar_reference` (same >12-bit
  float-accumulating vertical path as above, SSE side)
- `neon::plane_u16::neon_row_matches_scalar_reference` /
  `neon_rows_4_matches_scalar_reference`
- `neon::vertical_s16_hb::neon_vertical_hb_matches_scalar_reference`
  (signed 16-bit high-bit-depth vertical path)
- `neon::vertical_u16_hb::neon_vertical_hb_matches_scalar_reference`

Fix: audit the rounding/shift step in each affected kernel's final
narrow-and-round sequence and confirm it matches the scalar reference's
round-half-up-from-a-fixed-bias behavior bit-for-bit, paying particular
attention to the `f32`-accumulating "high bit depth" (>12-bit) vertical paths
and the AVX2 FMA-vs-default split, since those are where the divergence
concentrates.

## NEON RDM (rounding doubling multiply) rounding divergence

Same symptom class as above (small ±1 divergence, `Bilinear`-reproducible
here), specific to the RDM-accelerated NEON paths:

- `neon::plane_u8_rdm::neon_rdm_row_matches_scalar_reference` /
  `neon_rdm_rows_4_matches_scalar_reference`
- `neon::cbcr8_rdm::neon_rdm_row_matches_scalar_reference` /
  `neon_rdm_rows_4_matches_scalar_reference`

Fix: check the `vqrdmlahq_*` (saturating rounding doubling multiply-accumulate)
usage in these kernels for a rounding-bias or lane-fold difference vs. the
scalar Q15 reference.

## NEON `rgba_f16` horizontal convolution - confirmed real crash, not just a test failure

- **`convolve_horizontal_rgba_neon_row_one_f16`** (`src/neon/rgba_f16.rs:215`):
  the tail loop reads `while jx <= bounds.size` instead of `while jx <
  bounds.size` (every sibling loop in the same file, including the 4-row
  variant, correctly uses `<`). This performs one extra iteration, reading a
  source pixel and a filter weight **one past the valid range**.
  - **This is not a hypothetical bug.** Running the comparison test under
    `qemu-aarch64` with debug assertions enabled produces an actual
    `unsafe precondition(s) violated: slice::get_unchecked requires that the
    index is within the slice` panic that **aborts the whole process**
    (`SIGABRT`), confirming real out-of-bounds memory access, not just a
    wrong numeric result.
  - Test: `neon::rgba_f16::neon_row_matches_scalar_reference` (marked
    `#[ignore]` specifically because leaving it enabled crashes the entire
    test binary instead of just failing one test).
  - Fix: change `jx <= bounds.size` to `jx < bounds.size` at
    `src/neon/rgba_f16.rs:215` (one-line fix, mechanical - this is the
    highest-priority item in this list since it's a real memory-safety bug,
    not just an output-precision difference).

## SVE2

- **`sve_convolve_horizontal_rgb_neon_rows_4_dot`**
  (`src/sve2/rgb_u8_dot.rs`): the single-row variant of this function is
  bit-exact, but the 4-row variant produces a **degenerate, repeating
  pattern** for rows 1-3 of every 4-row batch - each output pixel comes out
  with `R == G == B`, and the same short sequence of values repeats three
  times across what should be three different rows of real image data.
  - Likely cause: `svext_u8::<4>`, `svext_u8::<8>`, `svext_u8::<12>`
    (`src/sve2/rgb_u8_dot.rs:170-172`) extract rows 1/2/3 from the packed
    4-row result using **hardcoded byte offsets sized for a 128-bit SVE
    vector length**. SVE/SVE2 is explicitly a *scalable* vector ISA - the
    real vector length is queried at runtime via `svcntb()`/`svcnth()` and
    can be wider than 128 bits on real hardware (and apparently is, by
    default, under this project's QEMU setup, given the test fails). When
    the true vector length isn't 128 bits, fixed `<4>/<8>/<12>` byte
    extents no longer land on row boundaries, corrupting rows 1-3.
  - This matches a *previously flagged, purely static-analysis* finding
    about this exact function (hardcoded vector-length assumptions in SVE2
    code) - this test run is the first **empirical confirmation** that it's
    a real, live bug, not just a theoretical code-smell.
  - Test: `sve2::rgb_u8_dot::sve2_rows_4_matches_scalar_reference`.
  - Fix: query the actual vector length (`svcntb()`) instead of hardcoding
    `4`/`8`/`12`, and rework the row-extraction logic to be correct for any
    SVE2 implementation width, not just 128-bit.

- **`convolve_vertical_sve2_u16_dot`** (`src/sve2/vertical_u16_dot.rs`): same
  small ±1 rounding divergence class as the NEON/AVX2/SSE u16 findings above.
  - Test: `sve2::vertical_u16_dot::sve2_vertical_matches_scalar_reference`.
  - Fix: same investigation as the "small (±1) fixed-point rounding
    divergences" section above - check whether this shares the same root
    cause.

Miri

Miri: memory-safety UB (found by cross-interpreting NEON with Miri, no ARM hardware or QEMU needed)

Running the same backend_comparison_tests under Miri (cargo +nightly miri test --target aarch64-unknown-linux-gnu ..., see just miri-neon) surfaces
real, language-level Undefined Behavior that QEMU can't detect, because
QEMU just executes real instructions on real (or emulated) hardware, which
tends to tolerate things the Rust abstract machine does not (over-aligned
allocators, forgiving unaligned-load handling on modern ARM cores, etc.).
Miri instead checks against the actual memory model, so these findings are
worth taking seriously even though nothing has been observed to crash on
real hardware yet.

Important caveat: Miri does not implement every NEON intrinsic this crate
uses. When it hits one it doesn't support, the WHOLE test binary aborts with
error: unsupported operation: can't call LLVM intrinsic ... - this is a
Miri limitation, not a bug, and it means that module simply cannot be
checked this way at all right now. Do not confuse this with
error: Undefined Behavior: ..., which always means a real, actionable bug.
just miri-neon runs module-by-module specifically so one module's abort
(for either reason) doesn't hide every module after it.

A 15-module sample (not exhaustive - there are ~70 files under src/neon/)
came back:

Module	Result
alpha_f32, rgba_f32	clean
plane_u16	clean (only the already-known rounding bug, ignored)
rgb_f32	real UB - Stacked Borrows violation
vertical_f32	real UB - out-of-bounds read
rgba_u8, rgb_u8, plane_u8, cbcr8, rgb_u16	real UB - misaligned pointer cast (same root cause, 5 modules)
plane_f32, rgba_u16, alpha_u16, vertical_u8, vertical_u16	Miri tooling gap (faddv/fcvtau/urshl/sqshrun unsupported) - not checked

Misaligned pointer cast in shared tail-pixel-load helpers (5+ modules)

src/neon/utils.rs has four small helpers used across the u8 NEON kernels to
load a partial (1-3 byte) tail pixel by reinterpreting a *const u8 as a
wider pointer type and using a single-lane NEON load intrinsic on it:

// src/neon/utils.rs:163-166
pub(crate) unsafe fn load_4b_as_u16x4(src_ptr: *const u8) -> uint16x4_t {
    unsafe {
        let j = vreinterpret_u8_u32(vld1_lane_u32::<0>(src_ptr as *const u32, vdup_n_u32(0)));
        ...

(load_3b_as_u16x4, load_3b_as_u8x16, load_4b_as_u8x8 at lines 145, 155,
172 all do the same thing at u16/u32 granularity.) vld1_lane_u32
requires its source pointer to be 4-byte aligned; a *const u8 computed from
an arbitrary byte offset into a &[u8] buffer offers no such guarantee, so
this is Undefined Behavior regardless of whether it happens to work on a
given allocator/CPU. Confirmed triggering it: rgba_u8, rgb_u8, plane_u8
(alignment 4 required), cbcr8 (alignment 2 required), rgb_u16 (alignment
4 required) - likely also rgb_u8_dot.rs, rgb_u8_sqrdml.rs,
rgba_u8_rdm.rs (confirmed callers of the same helpers, not yet run under
Miri individually).

• Fix: replace the aligned lane-load with an explicitly unaligned read, e.g.
  u32::from_ne_bytes(src_ptr.cast::<[u8; 4]>().read()) (or
  core::ptr::read_unaligned) to build the scalar value, then
  vdup_n_u32/vreinterpret_* as before - this produces the identical
  numeric result without requiring alignment the caller can't guarantee.

src/neon/rgb_f32.rs:76 - Stacked Borrows violation (reads past a narrowed reference's valid range)

// src/neon/rgb_f32.rs:76
let rgb_pixel_1 = unsafe { vld1q_f32(src_ptr.get_unchecked(2)) };

vld1q_f32 reads 4 f32s (16 bytes) starting at element offset 2, but the
reference obtained via get_unchecked(2) only carries provenance for a
narrower range - Miri traces this to an invalid access at the tail of a
row. Same family of bug as the confirmed QEMU crash in
src/neon/rgba_f16.rs:215 (an accelerated kernel reading a fixed SIMD width
past what the current tap/column actually has remaining) - likely the same
underlying "assumes a full vector's worth of data is always available"
mistake, here for f32 RGB rather than f16 RGBA.

• Fix: bound this load to the actual remaining width before issuing it (mirror
  whatever tail-handling the working sibling kernels in the same file already
  do correctly for smaller widths).

src/neon/vertical_f32.rs - genuine out-of-bounds read for narrow row widths

error: Undefined Behavior: memory access failed: attempting to access 16 bytes,
but got alloc+0x400 which is at or beyond the end of the allocation of size 1024 bytes
  --> .../vld1q_f32
      neon::utils::xvld1q_f32_x2 (src/neon/utils.rs:83)
      neon::vertical_f32::convolve_vertical_part_neon_4_f32 (src/neon/vertical_f32.rs:268)

xvld1q_f32_x2 loads two consecutive float32x4_t (8 floats / 32 bytes) at
once. convolve_vertical_part_neon_4_f32 calls it unconditionally even when
the row is narrower than 8 columns (reproduced with row width 4), so for
the last source row of a narrow image the second vld1q_f32 inside the x2
helper reads past the end of the entire source buffer, not just past the
current row - a genuine out-of-bounds read, not a test-sizing artifact (the
test buffer is sized generously at src_stride * in_height, matching how
every other vertical test in this suite sizes its buffer).

• Fix: gate the 8-wide (x2) load path on the actual remaining column count,
  falling back to a single 4-wide (or narrower) load for rows that don't have
  8 columns left, the same way the file's own narrower tail loops already do
  for other widths.

[awxkee]

Most of found UBs "Misaligned pointer cast" are false-positive bugs that come from incorrect debug implementations in Rust itself. I reported them a while ago, but currently no solution rust-lang/stdarch#2198. In release codegen seems to be correct.

[awxkee]

There is one real bug which I fixed in '#162' and not "so good" thing with stacked borrows in it.

However, some things as alpha premultiplication need a bit more investigation.