wide/

f32x4_.rs

use super::*;

pick! {
  if #[cfg(target_feature="sse")] {
    #[derive(Default, Clone, Copy, PartialEq)]
    #[repr(C, align(16))]
    pub struct f32x4 { pub(crate) sse: m128 }
  } else if #[cfg(target_feature="simd128")] {
    use core::arch::wasm32::*;

    #[derive(Clone, Copy)]
    #[repr(transparent)]
    pub struct f32x4 { pub(crate) simd: v128 }

    impl Default for f32x4 {
      fn default() -> Self {
        Self::splat(0.0)
      }
    }

    impl PartialEq for f32x4 {
      fn eq(&self, other: &Self) -> bool {
        u32x4_all_true(f32x4_eq(self.simd, other.simd))
      }
    }
  } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
    use core::arch::aarch64::*;
    #[repr(C)]
    #[derive(Copy, Clone)]
    pub struct f32x4 { pub(crate) neon : float32x4_t }

    impl Default for f32x4 {
      #[inline]
      fn default() -> Self {
        unsafe { Self { neon: vdupq_n_f32(0.0)} }
      }
    }

    impl PartialEq for f32x4 {
      #[inline]
      fn eq(&self, other: &Self) -> bool {
        unsafe { vminvq_u32(vceqq_f32(self.neon, other.neon))==u32::MAX }
      }

    }
    } else {
    #[derive(Default, Clone, Copy, PartialEq)]
    #[repr(C, align(16))]
    pub struct f32x4 { pub(crate) arr: [f32;4] }
  }
}

macro_rules! const_f32_as_f32x4 {
  ($i:ident, $f:expr) => {
    #[allow(non_upper_case_globals)]
    pub const $i: f32x4 = f32x4::new([$f; 4]);
  };
}

impl f32x4 {
  const_f32_as_f32x4!(ONE, 1.0);
  const_f32_as_f32x4!(ZERO, 0.0);
  const_f32_as_f32x4!(HALF, 0.5);
  const_f32_as_f32x4!(EPSILON, f32::EPSILON);
  const_f32_as_f32x4!(MIN, f32::MIN);
  const_f32_as_f32x4!(MIN_POSITIVE, f32::MIN_POSITIVE);
  const_f32_as_f32x4!(MAX, f32::MAX);
  const_f32_as_f32x4!(NAN, f32::NAN);
  const_f32_as_f32x4!(INFINITY, f32::INFINITY);
  const_f32_as_f32x4!(NEG_INFINITY, f32::NEG_INFINITY);
  const_f32_as_f32x4!(E, core::f32::consts::E);
  const_f32_as_f32x4!(FRAC_1_PI, core::f32::consts::FRAC_1_PI);
  const_f32_as_f32x4!(FRAC_2_PI, core::f32::consts::FRAC_2_PI);
  const_f32_as_f32x4!(FRAC_2_SQRT_PI, core::f32::consts::FRAC_2_SQRT_PI);
  const_f32_as_f32x4!(FRAC_1_SQRT_2, core::f32::consts::FRAC_1_SQRT_2);
  const_f32_as_f32x4!(FRAC_PI_2, core::f32::consts::FRAC_PI_2);
  const_f32_as_f32x4!(FRAC_PI_3, core::f32::consts::FRAC_PI_3);
  const_f32_as_f32x4!(FRAC_PI_4, core::f32::consts::FRAC_PI_4);
  const_f32_as_f32x4!(FRAC_PI_6, core::f32::consts::FRAC_PI_6);
  const_f32_as_f32x4!(FRAC_PI_8, core::f32::consts::FRAC_PI_8);
  const_f32_as_f32x4!(LN_2, core::f32::consts::LN_2);
  const_f32_as_f32x4!(LN_10, core::f32::consts::LN_10);
  const_f32_as_f32x4!(LOG2_E, core::f32::consts::LOG2_E);
  const_f32_as_f32x4!(LOG10_E, core::f32::consts::LOG10_E);
  const_f32_as_f32x4!(LOG10_2, core::f32::consts::LOG10_2);
  const_f32_as_f32x4!(LOG2_10, core::f32::consts::LOG2_10);
  const_f32_as_f32x4!(PI, core::f32::consts::PI);
  const_f32_as_f32x4!(SQRT_2, core::f32::consts::SQRT_2);
  const_f32_as_f32x4!(TAU, core::f32::consts::TAU);
}

unsafe impl Zeroable for f32x4 {}
unsafe impl Pod for f32x4 {}

impl AlignTo for f32x4 {
  type Elem = f32;
}

impl Add for f32x4 {
  type Output = Self;
  #[inline]
  fn add(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: add_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_add(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe { Self { neon: vaddq_f32(self.neon, rhs.neon) } }
      } else {
        Self { arr: [
          self.arr[0] + rhs.arr[0],
          self.arr[1] + rhs.arr[1],
          self.arr[2] + rhs.arr[2],
          self.arr[3] + rhs.arr[3],
        ]}
      }
    }
  }
}

impl Sub for f32x4 {
  type Output = Self;
  #[inline]
  fn sub(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: sub_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_sub(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vsubq_f32(self.neon, rhs.neon) }}
      } else {
        Self { arr: [
          self.arr[0] - rhs.arr[0],
          self.arr[1] - rhs.arr[1],
          self.arr[2] - rhs.arr[2],
          self.arr[3] - rhs.arr[3],
        ]}
      }
    }
  }
}

impl Mul for f32x4 {
  type Output = Self;
  #[inline]
  fn mul(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: mul_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_mul(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vmulq_f32(self.neon, rhs.neon) }}
      } else {
        Self { arr: [
          self.arr[0] * rhs.arr[0],
          self.arr[1] * rhs.arr[1],
          self.arr[2] * rhs.arr[2],
          self.arr[3] * rhs.arr[3],
        ]}
      }
    }
  }
}

impl Div for f32x4 {
  type Output = Self;
  #[inline]
  fn div(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: div_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_div(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vdivq_f32(self.neon, rhs.neon) }}
      } else {
        Self { arr: [
          self.arr[0] / rhs.arr[0],
          self.arr[1] / rhs.arr[1],
          self.arr[2] / rhs.arr[2],
          self.arr[3] / rhs.arr[3],
        ]}
      }
    }
  }
}

impl Rem for f32x4 {
  type Output = Self;
  #[inline]
  fn rem(self, rhs: Self) -> Self::Output {
    Self::new([
      self.to_array()[0] % rhs.to_array()[0],
      self.to_array()[1] % rhs.to_array()[1],
      self.to_array()[2] % rhs.to_array()[2],
      self.to_array()[3] % rhs.to_array()[3],
    ])
  }
}

impl Neg for f32x4 {
  type Output = Self;
  #[inline]
  fn neg(self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: bitxor_m128(self.sse, Self::splat(-0.0).sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_neg(self.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vnegq_f32(self.neon) }}
      } else {
        Self { arr: [
          -self.arr[0],
          -self.arr[1],
          -self.arr[2],
          -self.arr[3],
        ]}
      }
    }
  }
}

impl Add<f32> for f32x4 {
  type Output = Self;
  #[inline]
  fn add(self, rhs: f32) -> Self::Output {
    self.add(Self::splat(rhs))
  }
}

impl Sub<f32> for f32x4 {
  type Output = Self;
  #[inline]
  fn sub(self, rhs: f32) -> Self::Output {
    self.sub(Self::splat(rhs))
  }
}

impl Mul<f32> for f32x4 {
  type Output = Self;
  #[inline]
  fn mul(self, rhs: f32) -> Self::Output {
    self.mul(Self::splat(rhs))
  }
}

impl Div<f32> for f32x4 {
  type Output = Self;
  #[inline]
  fn div(self, rhs: f32) -> Self::Output {
    self.div(Self::splat(rhs))
  }
}

impl Rem<f32> for f32x4 {
  type Output = Self;
  #[inline]
  fn rem(self, rhs: f32) -> Self::Output {
    self.rem(Self::splat(rhs))
  }
}

impl Add<f32x4> for f32 {
  type Output = f32x4;
  #[inline]
  fn add(self, rhs: f32x4) -> Self::Output {
    f32x4::splat(self).add(rhs)
  }
}

impl Sub<f32x4> for f32 {
  type Output = f32x4;
  #[inline]
  fn sub(self, rhs: f32x4) -> Self::Output {
    f32x4::splat(self).sub(rhs)
  }
}

impl Mul<f32x4> for f32 {
  type Output = f32x4;
  #[inline]
  fn mul(self, rhs: f32x4) -> Self::Output {
    f32x4::splat(self).mul(rhs)
  }
}

impl Div<f32x4> for f32 {
  type Output = f32x4;
  #[inline]
  fn div(self, rhs: f32x4) -> Self::Output {
    f32x4::splat(self).div(rhs)
  }
}

impl Rem<f32x4> for f32 {
  type Output = f32x4;
  #[inline]
  fn rem(self, rhs: f32x4) -> Self::Output {
    f32x4::splat(self).rem(rhs)
  }
}

impl BitAnd for f32x4 {
  type Output = Self;
  #[inline]
  fn bitand(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: bitand_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: v128_and(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vandq_u32(vreinterpretq_u32_f32(self.neon), vreinterpretq_u32_f32(rhs.neon))) }}
      } else {
        Self { arr: [
          f32::from_bits(self.arr[0].to_bits() & rhs.arr[0].to_bits()),
          f32::from_bits(self.arr[1].to_bits() & rhs.arr[1].to_bits()),
          f32::from_bits(self.arr[2].to_bits() & rhs.arr[2].to_bits()),
          f32::from_bits(self.arr[3].to_bits() & rhs.arr[3].to_bits()),
        ]}
      }
    }
  }
}

impl BitOr for f32x4 {
  type Output = Self;
  #[inline]
  fn bitor(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: bitor_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: v128_or(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vorrq_u32(vreinterpretq_u32_f32(self.neon), vreinterpretq_u32_f32(rhs.neon))) }}
      } else {
        Self { arr: [
          f32::from_bits(self.arr[0].to_bits() | rhs.arr[0].to_bits()),
          f32::from_bits(self.arr[1].to_bits() | rhs.arr[1].to_bits()),
          f32::from_bits(self.arr[2].to_bits() | rhs.arr[2].to_bits()),
          f32::from_bits(self.arr[3].to_bits() | rhs.arr[3].to_bits()),
        ]}
      }
    }
  }
}

impl BitXor for f32x4 {
  type Output = Self;
  #[inline]
  fn bitxor(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: bitxor_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: v128_xor(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(veorq_u32(vreinterpretq_u32_f32(self.neon), vreinterpretq_u32_f32(rhs.neon))) }}
      } else {
        Self { arr: [
          f32::from_bits(self.arr[0].to_bits() ^ rhs.arr[0].to_bits()),
          f32::from_bits(self.arr[1].to_bits() ^ rhs.arr[1].to_bits()),
          f32::from_bits(self.arr[2].to_bits() ^ rhs.arr[2].to_bits()),
          f32::from_bits(self.arr[3].to_bits() ^ rhs.arr[3].to_bits()),
        ]}
      }
    }
  }
}

#[expect(deprecated)]
impl CmpEq for f32x4 {
  type Output = Self;
  #[inline]
  fn simd_eq(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: cmp_eq_mask_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_eq(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vceqq_f32(self.neon, rhs.neon)) }}
      } else {
        Self { arr: [
          if self.arr[0] == rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[1] == rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[2] == rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[3] == rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
        ]}
      }
    }
  }
}

#[expect(deprecated)]
impl CmpGe for f32x4 {
  type Output = Self;
  #[inline]
  fn simd_ge(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: cmp_ge_mask_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_ge(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vcgeq_f32(self.neon, rhs.neon)) }}
      } else {
        Self { arr: [
          if self.arr[0] >= rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[1] >= rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[2] >= rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[3] >= rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
        ]}
      }
    }
  }
}

#[expect(deprecated)]
impl CmpGt for f32x4 {
  type Output = Self;
  #[inline]
  fn simd_gt(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: cmp_gt_mask_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_gt(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vcgtq_f32(self.neon, rhs.neon)) }}
      } else {
        Self { arr: [
          if self.arr[0] > rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[1] > rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[2] > rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[3] > rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
        ]}
      }
    }
  }
}

#[expect(deprecated)]
impl CmpNe for f32x4 {
  type Output = Self;
  #[inline]
  fn simd_ne(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: cmp_neq_mask_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_ne(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vmvnq_u32(vceqq_f32(self.neon, rhs.neon))) }}
      } else {
        Self { arr: [
          if self.arr[0] != rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[1] != rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[2] != rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[3] != rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
        ]}
      }
    }
  }
}

#[expect(deprecated)]
impl CmpLe for f32x4 {
  type Output = Self;
  #[inline]
  fn simd_le(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: cmp_le_mask_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_le(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vcleq_f32(self.neon, rhs.neon)) }}
      } else {
        Self { arr: [
          if self.arr[0] <= rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[1] <= rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[2] <= rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[3] <= rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
        ]}
      }
    }
  }
}

#[expect(deprecated)]
impl CmpLt for f32x4 {
  type Output = Self;
  #[inline]
  fn simd_lt(self, rhs: Self) -> Self::Output {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: cmp_lt_mask_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_lt(self.simd, rhs.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vcltq_f32(self.neon, rhs.neon)) }}
      } else {
        Self { arr: [
          if self.arr[0] < rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[1] < rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[2] < rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[3] < rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
        ]}
      }
    }
  }
}

impl f32x4 {
  #[inline]
  #[must_use]
  pub const fn new(array: [f32; 4]) -> Self {
    unsafe { core::mem::transmute(array) }
  }

  simd_comparison_fns!();

  #[inline]
  #[must_use]
  pub fn blend(self, t: Self, f: Self) -> Self {
    pick! {
      if #[cfg(target_feature="sse4.1")] {
        Self { sse: blend_varying_m128(f.sse, t.sse, self.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: v128_bitselect(t.simd, f.simd, self.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vbslq_f32(vreinterpretq_u32_f32(self.neon), t.neon, f.neon) }}
      } else {
        generic_bit_blend(self, t, f)
      }
    }
  }
  #[inline]
  #[must_use]
  pub fn abs(self) -> Self {
    pick! {
      if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_abs(self.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vabsq_f32(self.neon) }}
      } else {
        let non_sign_bits = f32x4::from(f32::from_bits(i32::MAX as u32));
        self & non_sign_bits
      }
    }
  }

  #[inline]
  #[must_use]
  pub fn signum(self) -> Self {
    let result = Self::ONE | self & -Self::ZERO;

    self.is_nan().blend(self, result)
  }

  #[inline]
  #[must_use]
  pub fn floor(self) -> Self {
    pick! {
      if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_floor(self.simd) }
      } else if #[cfg(target_feature="sse4.1")] {
        Self { sse: floor_m128(self.sse) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vrndmq_f32(self.neon) }}
      } else if #[cfg(feature="std")] {
        let base: [f32; 4] = cast(self);
        cast(base.map(|val| val.floor()))
      } else {
        let base: [f32; 4] = cast(self);
        let rounded: [f32; 4] = cast(self.round());
        cast([
          if base[0] < rounded[0] { rounded[0] - 1.0 } else { rounded[0] },
          if base[1] < rounded[1] { rounded[1] - 1.0 } else { rounded[1] },
          if base[2] < rounded[2] { rounded[2] - 1.0 } else { rounded[2] },
          if base[3] < rounded[3] { rounded[3] - 1.0 } else { rounded[3] },
        ])
      }
    }
  }
  #[inline]
  #[must_use]
  pub fn ceil(self) -> Self {
    pick! {
      if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_ceil(self.simd) }
      } else if #[cfg(target_feature="sse4.1")] {
        Self { sse: ceil_m128(self.sse) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vrndpq_f32(self.neon) }}
      } else if #[cfg(feature="std")] {
        let base: [f32; 4] = cast(self);
        cast(base.map(|val| val.ceil()))
      } else {
        let base: [f32; 4] = cast(self);
        let rounded: [f32; 4] = cast(self.round());
        cast([
          if base[0] > rounded[0] { rounded[0] + 1.0 } else { rounded[0] },
          if base[1] > rounded[1] { rounded[1] + 1.0 } else { rounded[1] },
          if base[2] > rounded[2] { rounded[2] + 1.0 } else { rounded[2] },
          if base[3] > rounded[3] { rounded[3] + 1.0 } else { rounded[3] },
        ])
      }
    }
  }

  /// Calculates the lanewise maximum of both vectors. This is a faster
  /// implementation than `max`, but it doesn't specify any behavior if NaNs are
  /// involved.
  #[inline]
  #[must_use]
  pub fn fast_max(self, rhs: Self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: max_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self {
          simd: f32x4_pmax(self.simd, rhs.simd),
        }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vmaxq_f32(self.neon, rhs.neon) }}
      } else {
        Self { arr: [
          if self.arr[0] < rhs.arr[0] { rhs.arr[0] } else { self.arr[0] },
          if self.arr[1] < rhs.arr[1] { rhs.arr[1] } else { self.arr[1] },
          if self.arr[2] < rhs.arr[2] { rhs.arr[2] } else { self.arr[2] },
          if self.arr[3] < rhs.arr[3] { rhs.arr[3] } else { self.arr[3] },
        ]}
      }
    }
  }

  /// Calculates the lanewise maximum of both vectors. If either lane is NaN,
  /// the other lane gets chosen. Use `fast_max` for a faster implementation
  /// that doesn't handle NaNs.
  #[inline]
  #[must_use]
  pub fn max(self, rhs: Self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        // max_m128 seems to do rhs < self ? self : rhs. So if there's any NaN
        // involved, it chooses rhs, so we need to specifically check rhs for
        // NaN.
        rhs.is_nan().blend(self, Self { sse: max_m128(self.sse, rhs.sse) })
      } else if #[cfg(target_feature="simd128")] {
        // WASM has two max intrinsics:
        // - max: This propagates NaN, that's the opposite of what we need.
        // - pmax: This is defined as self < rhs ? rhs : self, which basically
        //   chooses self if either is NaN.
        //
        // pmax is what we want, but we need to specifically check self for NaN.
        Self {
          simd: v128_bitselect(
            rhs.simd,
            f32x4_pmax(self.simd, rhs.simd),
            f32x4_ne(self.simd, self.simd), // NaN check
          )
        }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vmaxnmq_f32(self.neon, rhs.neon) }}
      } else {
        Self { arr: [
          self.arr[0].max(rhs.arr[0]),
          self.arr[1].max(rhs.arr[1]),
          self.arr[2].max(rhs.arr[2]),
          self.arr[3].max(rhs.arr[3]),
        ]}
      }
    }
  }

  /// Calculates the lanewise minimum of both vectors. This is a faster
  /// implementation than `min`, but it doesn't specify any behavior if NaNs are
  /// involved.
  #[inline]
  #[must_use]
  pub fn fast_min(self, rhs: Self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: min_m128(self.sse, rhs.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self {
          simd: f32x4_pmin(self.simd, rhs.simd),
        }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vminq_f32(self.neon, rhs.neon) }}
      } else {
        Self { arr: [
          if self.arr[0] < rhs.arr[0] { self.arr[0] } else { rhs.arr[0] },
          if self.arr[1] < rhs.arr[1] { self.arr[1] } else { rhs.arr[1] },
          if self.arr[2] < rhs.arr[2] { self.arr[2] } else { rhs.arr[2] },
          if self.arr[3] < rhs.arr[3] { self.arr[3] } else { rhs.arr[3] },
        ]}
      }
    }
  }

  /// Calculates the lanewise minimum of both vectors. If either lane is NaN,
  /// the other lane gets chosen. Use `fast_min` for a faster implementation
  /// that doesn't handle NaNs.
  #[inline]
  #[must_use]
  pub fn min(self, rhs: Self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        // min_m128 seems to do self < rhs ? self : rhs. So if there's any NaN
        // involved, it chooses rhs, so we need to specifically check rhs for
        // NaN.
        rhs.is_nan().blend(self, Self { sse: min_m128(self.sse, rhs.sse) })
      } else if #[cfg(target_feature="simd128")] {
        // WASM has two min intrinsics:
        // - min: This propagates NaN, that's the opposite of what we need.
        // - pmin: This is defined as rhs < self ? rhs : self, which basically
        //   chooses self if either is NaN.
        //
        // pmin is what we want, but we need to specifically check self for NaN.
        Self {
          simd: v128_bitselect(
            rhs.simd,
            f32x4_pmin(self.simd, rhs.simd),
            f32x4_ne(self.simd, self.simd), // NaN check
          )
        }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vminnmq_f32(self.neon, rhs.neon) }}
      } else {
        Self { arr: [
          self.arr[0].min(rhs.arr[0]),
          self.arr[1].min(rhs.arr[1]),
          self.arr[2].min(rhs.arr[2]),
          self.arr[3].min(rhs.arr[3]),
        ]}
      }
    }
  }

  /// Restrict a value to a certain interval unless it is NaN.
  ///
  /// If `self` is NaN, or `min` is NaN, or `max` is NaN, the result is NaN.
  /// If `min > max`, the result is `min`, since `fast_max(min)` dominates.
  #[inline]
  #[must_use]
  pub fn clamp(self, min: Self, max: Self) -> Self {
    let is_nan = self.is_nan() | min.is_nan() | max.is_nan();
    let clamped = self.fast_min(max).fast_max(min);
    is_nan.blend(Self::splat(f32::NAN), clamped)
  }

  /// Restrict a value to a certain interval unless it is NaN.
  ///
  /// Avoids NaN detection; same speed as the old `clamp` prior to IEEE 754-2019
  /// compliance. Does not specify any
  /// behavior if NaNs are involved, and if `min > max` the result is
  /// unspecified.
  #[inline]
  #[must_use]
  pub fn fast_clamp(self, min: Self, max: Self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        // For both `min_m128` and `max_m128` if any input is NaN, `rhs` gets
        // chosen. For `self` to be chosen, `self` must be the second argument.
        Self { sse: min_m128(max.sse, max_m128(min.sse, self.sse)) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_min(f32x4_max(self.simd, min.simd), max.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
        unsafe { Self { neon: vminq_f32(vmaxq_f32(self.neon, min.neon), max.neon) } }
      } else {
        // The standard library does not have NaN propagating `min` and `max`
        // functions.
        let mut result = self;
        result = result.simd_lt(min).blend(min, self);
        result = result.simd_gt(max).blend(max, self);
        result
      }
    }
  }

  #[inline]
  #[must_use]
  pub fn midpoint(self, other: Self) -> Self {
    (self + other) * 0.5
  }

  #[inline]
  #[must_use]
  pub fn is_nan(self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: cmp_unord_mask_m128(self.sse, self.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_ne(self.simd, self.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vreinterpretq_f32_u32(vmvnq_u32(vceqq_f32(self.neon, self.neon))) }}
      } else {
        Self { arr: [
          if self.arr[0].is_nan() { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[1].is_nan() { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[2].is_nan() { f32::from_bits(u32::MAX) } else { 0.0 },
          if self.arr[3].is_nan() { f32::from_bits(u32::MAX) } else { 0.0 },
        ]}
      }
    }
  }
  #[inline]
  #[must_use]
  pub fn is_finite(self) -> Self {
    let shifted_exp_mask = u32x4::from(0xFF000000);
    let u: u32x4 = cast(self);
    let shift_u = u << 1_u64;
    let out = !(shift_u & shifted_exp_mask).simd_eq(shifted_exp_mask);
    cast(out)
  }
  #[inline]
  #[must_use]
  pub fn is_inf(self) -> Self {
    let shifted_inf = u32x4::from(0xFF000000);
    let u: u32x4 = cast(self);
    let shift_u = u << 1_u64;
    let out = (shift_u).simd_eq(shifted_inf);
    cast(out)
  }

  #[inline]
  #[must_use]
  pub fn round(self) -> Self {
    pick! {
      if #[cfg(target_feature="sse4.1")] {
        Self { sse: round_m128::<{round_op!(Nearest)}>(self.sse) }
      } else if #[cfg(target_feature="sse2")] {
        let mi: m128i = convert_to_i32_m128i_from_m128(self.sse);
        let f: f32x4 = f32x4 { sse: convert_to_m128_from_i32_m128i(mi) };
        let i: i32x4 = cast(mi);
        let mask: f32x4 = cast(i.simd_eq(i32x4::from(0x80000000_u32 as i32)));
        mask.blend(self, f)
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_nearest(self.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vrndnq_f32(self.neon) }}
      } else {
        // Note(Lokathor): This software fallback is probably very slow compared
        // to having a hardware option available, even just the sse2 version is
        // better than this. Oh well.
        let to_int = f32x4::from(1.0 / f32::EPSILON);
        let u: u32x4 = cast(self);
        let e: i32x4 = cast((u >> 23) & u32x4::from(0xff));
        let mut y: f32x4;

        let no_op_magic = i32x4::from(0x7f + 23);
        let no_op_mask: f32x4 = cast(e.simd_gt(no_op_magic) | e.simd_eq(no_op_magic));
        let no_op_val: f32x4 = self;

        let zero_magic = i32x4::from(0x7f - 1);
        let zero_mask: f32x4 = cast(e.simd_lt(zero_magic));
        let zero_val: f32x4 = self * f32x4::from(0.0);

        let neg_bit: f32x4 = cast(cast::<u32x4, i32x4>(u).simd_lt(i32x4::default()));
        let x: f32x4 = neg_bit.blend(-self, self);
        y = x + to_int - to_int - x;
        y = y.simd_gt(f32x4::from(0.5)).blend(
          y + x - f32x4::from(-1.0),
          y.simd_lt(f32x4::from(-0.5)).blend(y + x + f32x4::from(1.0), y + x),
        );
        y = neg_bit.blend(-y, y);

        no_op_mask.blend(no_op_val, zero_mask.blend(zero_val, y))
      }
    }
  }

  /// Rounds each lane into an integer. This is a faster implementation than
  /// `round_int`, but it doesn't handle out of range values or NaNs. For those
  /// values you get implementation defined behavior.
  #[inline]
  #[must_use]
  pub fn fast_round_int(self) -> i32x4 {
    pick! {
      if #[cfg(target_feature="sse2")] {
        cast(convert_to_i32_m128i_from_m128(self.sse))
      } else {
        self.round_int()
      }
    }
  }

  /// Rounds each lane into an integer. This saturates out of range values and
  /// turns NaNs into 0. Use `fast_round_int` for a faster implementation that
  /// doesn't handle out of range values or NaNs.
  #[inline]
  #[must_use]
  pub fn round_int(self) -> i32x4 {
    pick! {
      if #[cfg(target_feature="sse2")] {
        // Based on: https://github.com/v8/v8/blob/210987a552a2bf2a854b0baa9588a5959ff3979d/src/codegen/shared-ia32-x64/macro-assembler-shared-ia32-x64.h#L489-L504
        let non_nan_mask = self.simd_eq(self);
        let non_nan = self & non_nan_mask;
        let flip_to_max: i32x4 = cast(self.simd_ge(Self::splat(2147483648.0)));
        let cast: i32x4 = cast(convert_to_i32_m128i_from_m128(non_nan.sse));
        flip_to_max ^ cast
      } else if #[cfg(target_feature="simd128")] {
        cast(Self { simd: i32x4_trunc_sat_f32x4(f32x4_nearest(self.simd)) })
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        cast(unsafe {Self { neon: vreinterpretq_f32_s32(vcvtnq_s32_f32(self.neon)) }})
      } else {
        let rounded: [f32; 4] = cast(self.round());
        cast([
          rounded[0] as i32,
          rounded[1] as i32,
          rounded[2] as i32,
          rounded[3] as i32,
        ])
      }
    }
  }

  #[inline]
  #[must_use]
  pub fn trunc(self) -> Self {
    pick! {
      if #[cfg(target_feature="sse4.1")] {
        Self { sse: round_m128::<{round_op!(Zero)}>(self.sse) }
      } else if #[cfg(target_feature="sse2")] {
        // Ported from https://docs.rs/glam/latest/glam/f32/struct.Vec4.html#method.trunc
        // Based on https://github.com/microsoft/DirectXMath `XMVectorTruncate`
        let result: Self = cast(convert_to_m128_from_i32_m128i(truncate_m128_to_m128i(self.sse)));

        // Out of range values are either already round, infinite or NaN.
        let bounds_mask: Self = cast(cmp_lt_mask_i32_m128i(
            cast(self.abs()),
            set_splat_i32_m128i(8388608_f32.to_bits() as i32),
        ));

        // Reset the sign bit of the mask to preverse the sign of `self`.
        bounds_mask.abs().blend(result, self)
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_trunc(self.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe { Self { neon: vrndq_f32(self.neon) } }
      } else {
        let array: [f32; 4] = cast(self);
        let result: Self = cast([
          array[0] as i32 as f32,
          array[1] as i32 as f32,
          array[2] as i32 as f32,
          array[3] as i32 as f32,
        ]);

        // Out of range values are either already round, infinite or NaN.
        const BOUNDS_LIMIT: i32 = 8388608_f32.to_bits() as i32;
        let bounds_mask: Self = cast(cast::<f32x4, i32x4>(self.abs()).simd_lt(i32x4::splat(BOUNDS_LIMIT)));

        // Reset the sign bit of the mask to preverse the sign of `self`.
        bounds_mask.abs().blend(result, self)
      }
    }
  }

  /// Truncates each lane into an integer. This is a faster implementation than
  /// `trunc_int`, but it doesn't handle out of range values or NaNs. For those
  /// values you get implementation defined behavior.
  #[inline]
  #[must_use]
  pub fn fast_trunc_int(self) -> i32x4 {
    pick! {
      if #[cfg(target_feature="sse2")] {
        cast(truncate_m128_to_m128i(self.sse))
      } else {
        self.trunc_int()
      }
    }
  }

  /// Truncates each lane into an integer. This saturates out of range values
  /// and turns NaNs into 0. Use `fast_trunc_int` for a faster implementation
  /// that doesn't handle out of range values or NaNs.
  #[inline]
  #[must_use]
  pub fn trunc_int(self) -> i32x4 {
    pick! {
      if #[cfg(target_feature="sse2")] {
        // Based on: https://github.com/v8/v8/blob/210987a552a2bf2a854b0baa9588a5959ff3979d/src/codegen/shared-ia32-x64/macro-assembler-shared-ia32-x64.h#L489-L504
        let non_nan_mask = self.simd_eq(self);
        let non_nan = self & non_nan_mask;
        let flip_to_max: i32x4 = cast(self.simd_ge(Self::splat(2147483648.0)));
        let cast: i32x4 = cast(truncate_m128_to_m128i(non_nan.sse));
        flip_to_max ^ cast
      } else if #[cfg(target_feature="simd128")] {
        cast(Self { simd: i32x4_trunc_sat_f32x4(self.simd) })
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        cast(unsafe {Self { neon: vreinterpretq_f32_s32(vcvtq_s32_f32(self.neon)) }})
      } else {
        let n: [f32;4] = cast(self);
        cast([
          n[0] as i32,
          n[1] as i32,
          n[2] as i32,
          n[3] as i32,
        ])
      }
    }
  }

  #[inline]
  #[must_use]
  pub fn fract(self) -> Self {
    self - self.trunc()
  }

  /// Performs a multiply-add operation: `self * m + a`
  ///
  /// When hardware FMA support is available, this computes the result with a
  /// single rounding operation. Without FMA support, it falls back to separate
  /// multiply and add operations with two roundings.
  ///
  /// # Platform-specific behavior
  /// - On `x86`/`x86_64` with FMA: Uses `vfmadd` (single rounding, best
  ///   accuracy)
  /// - On ARM64 with NEON: Uses `vfmaq_f32` (single rounding, best accuracy)
  /// - Without FMA support: Uses `(self * m) + a` (two roundings)
  ///
  /// # Examples
  /// ```
  /// # use wide::f32x4;
  /// let a = f32x4::from([1.0, 2.0, 3.0, 4.0]);
  /// let b = f32x4::from([5.0, 6.0, 7.0, 8.0]);
  /// let c = f32x4::from([9.0, 10.0, 11.0, 12.0]);
  ///
  /// let result = a.mul_add(b, c);
  ///
  /// let expected = f32x4::from([14.0, 22.0, 32.0, 44.0]);
  /// assert_eq!(result, expected);
  /// ```
  #[inline]
  #[must_use]
  pub fn mul_add(self, m: Self, a: Self) -> Self {
    pick! {
      if #[cfg(all(target_feature="sse2",target_feature="fma"))] {
        Self { sse: fused_mul_add_m128(self.sse, m.sse, a.sse) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
        unsafe { Self { neon: vfmaq_f32(a.neon, self.neon, m.neon) } }
      } else {
        (self * m) + a
      }
    }
  }

  /// Performs a multiply-subtract operation: `self * m - s`
  ///
  /// When hardware FMA support is available, this computes the result with a
  /// single rounding operation. Without FMA support, it falls back to separate
  /// multiply and subtract operations with two roundings.
  ///
  /// # Platform-specific behavior
  /// - On `x86`/`x86_64` with FMA: Uses `vfmsub` (single rounding, best
  ///   accuracy)
  /// - On ARM64 with NEON: Uses `vfmaq_f32(-s, self, m)` (single rounding, best
  ///   accuracy)
  /// - Without FMA support: Uses `(self * m) - s` (two roundings)
  ///
  /// # Examples
  /// ```
  /// # use wide::f32x4;
  /// let a = f32x4::from([10.0, 20.0, 30.0, 40.0]);
  /// let b = f32x4::from([2.0, 3.0, 4.0, 5.0]);
  /// let c = f32x4::from([5.0, 10.0, 15.0, 20.0]);
  ///
  /// let result = a.mul_sub(b, c);
  ///
  /// let expected = f32x4::from([15.0, 50.0, 105.0, 180.0]);
  /// assert_eq!(result, expected);
  /// ```
  #[inline]
  #[must_use]
  pub fn mul_sub(self, m: Self, s: Self) -> Self {
    pick! {
      if #[cfg(all(target_feature="sse2",target_feature="fma"))] {
        Self { sse: fused_mul_sub_m128(self.sse, m.sse, s.sse) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
        unsafe { Self { neon: vfmaq_f32(vnegq_f32(s.neon), self.neon, m.neon) } }
      } else {
        (self * m) - s
      }
    }
  }

  /// Performs a negative multiply-add operation: `a - (self * m)`
  ///
  /// When hardware FMA support is available, this computes the result with a
  /// single rounding operation. Without FMA support, it falls back to separate
  /// operations with two roundings.
  ///
  /// # Platform-specific behavior
  /// - On `x86`/`x86_64` with FMA: Uses `vfnmadd` (single rounding, best
  ///   accuracy)
  /// - On ARM64 with NEON: Uses `vfmsq_f32` (single rounding, best accuracy)
  /// - Without FMA support: Uses `a - (self * m)` (two roundings)
  ///
  /// # Examples
  /// ```
  /// # use wide::f32x4;
  /// let a = f32x4::from([3.0, 4.0, 5.0, 6.0]);
  /// let b = f32x4::from([2.0, 2.0, 2.0, 2.0]);
  /// let c = f32x4::from([10.0, 20.0, 30.0, 40.0]);
  ///
  /// let result = a.mul_neg_add(b, c);
  ///
  /// let expected = f32x4::from([4.0, 12.0, 20.0, 28.0]);
  /// assert_eq!(result, expected);
  /// ```
  #[inline]
  #[must_use]
  pub fn mul_neg_add(self, m: Self, a: Self) -> Self {
    pick! {
      if #[cfg(all(target_feature="sse2",target_feature="fma"))] {
        Self { sse: fused_mul_neg_add_m128(self.sse, m.sse, a.sse) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
        unsafe { Self { neon: vfmsq_f32(a.neon, self.neon, m.neon) } }
      } else {
        a - (self * m)
      }
    }
  }

  /// Performs a negative multiply-subtract operation: `-(self * m) - s`
  ///
  /// When hardware FMA support is available, this computes the result with a
  /// single rounding operation. Without FMA support, it falls back to separate
  /// operations with two roundings.
  ///
  /// # Platform-specific behavior
  /// - On `x86`/`x86_64` with FMA: Uses `vfnmsub` (single rounding, best
  ///   accuracy)
  /// - On ARM64 with NEON: Uses `-(vfmaq_f32(s, self, m))` (single rounding,
  ///   best accuracy)
  /// - Without FMA support: Uses `-(self * m) - s` (two roundings)
  ///
  /// # Examples
  /// ```
  /// # use wide::f32x4;
  /// let a = f32x4::from([3.0, 4.0, 5.0, 6.0]);
  /// let b = f32x4::from([2.0, 2.0, 2.0, 2.0]);
  /// let c = f32x4::from([1.0, 2.0, 3.0, 4.0]);
  ///
  /// let result = a.mul_neg_sub(b, c);
  ///
  /// let expected = f32x4::from([-7.0, -10.0, -13.0, -16.0]);
  /// assert_eq!(result, expected);
  /// ```
  #[inline]
  #[must_use]
  pub fn mul_neg_sub(self, m: Self, s: Self) -> Self {
    pick! {
      if #[cfg(all(target_feature="sse2",target_feature="fma"))] {
        Self { sse: fused_mul_neg_sub_m128(self.sse, m.sse, s.sse) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
        unsafe { Self { neon: vnegq_f32(vfmaq_f32(s.neon, self.neon, m.neon)) } }
      } else {
        -(self * m) - s
      }
    }
  }

  #[inline]
  #[must_use]
  pub fn div_euclid(self, rhs: Self) -> Self {
    let q = (self / rhs).trunc();
    (self % rhs)
      .simd_lt(Self::ZERO)
      .blend(rhs.simd_gt(Self::ZERO).blend(q - Self::ONE, q + Self::ONE), q)
  }

  #[inline]
  #[must_use]
  pub fn rem_euclid(self, rhs: Self) -> Self {
    let r = self % rhs;
    r.simd_lt(Self::ZERO).blend(r + rhs.abs(), r)
  }

  #[inline]
  #[must_use]
  pub fn flip_signs(self, signs: Self) -> Self {
    self ^ (signs & Self::from(-0.0))
  }

  #[inline]
  #[must_use]
  pub fn copysign(self, sign: Self) -> Self {
    let magnitude_mask = Self::from(f32::from_bits(u32::MAX >> 1));
    (self & magnitude_mask) | (sign & Self::from(-0.0))
  }

  #[inline]
  pub fn asin_acos(self) -> (Self, Self) {
    // Based on the Agner Fog "vector class library":
    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
    const_f32_as_f32x4!(P4asinf, 4.2163199048E-2);
    const_f32_as_f32x4!(P3asinf, 2.4181311049E-2);
    const_f32_as_f32x4!(P2asinf, 4.5470025998E-2);
    const_f32_as_f32x4!(P1asinf, 7.4953002686E-2);
    const_f32_as_f32x4!(P0asinf, 1.6666752422E-1);

    let xa = self.abs();
    let big = xa.simd_ge(f32x4::splat(0.5));

    let x1 = f32x4::splat(0.5) * (f32x4::ONE - xa);
    let x2 = xa * xa;
    let x3 = big.blend(x1, x2);

    let xb = x1.sqrt();

    let x4 = big.blend(xb, xa);

    let z = polynomial_4!(x3, P0asinf, P1asinf, P2asinf, P3asinf, P4asinf);
    let z = z.mul_add(x3 * x4, x4);

    let z1 = z + z;

    // acos
    let z3 = self.simd_lt(f32x4::ZERO).blend(f32x4::PI - z1, z1);
    let z4 = f32x4::FRAC_PI_2 - z.flip_signs(self);
    let acos = big.blend(z3, z4);

    // asin
    let z3 = f32x4::FRAC_PI_2 - z1;
    let asin = big.blend(z3, z);
    let asin = asin.flip_signs(self);

    (asin, acos)
  }

  #[inline]
  pub fn asin(self) -> Self {
    // Based on the Agner Fog "vector class library":
    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
    const_f32_as_f32x4!(P4asinf, 4.2163199048E-2);
    const_f32_as_f32x4!(P3asinf, 2.4181311049E-2);
    const_f32_as_f32x4!(P2asinf, 4.5470025998E-2);
    const_f32_as_f32x4!(P1asinf, 7.4953002686E-2);
    const_f32_as_f32x4!(P0asinf, 1.6666752422E-1);

    let xa = self.abs();
    let big = xa.simd_ge(f32x4::splat(0.5));

    let x1 = f32x4::splat(0.5) * (f32x4::ONE - xa);
    let x2 = xa * xa;
    let x3 = big.blend(x1, x2);

    let xb = x1.sqrt();

    let x4 = big.blend(xb, xa);

    let z = polynomial_4!(x3, P0asinf, P1asinf, P2asinf, P3asinf, P4asinf);
    let z = z.mul_add(x3 * x4, x4);

    let z1 = z + z;

    // asin
    let z3 = f32x4::FRAC_PI_2 - z1;
    let asin = big.blend(z3, z);
    let asin = asin.flip_signs(self);

    asin
  }

  #[inline]
  #[must_use]
  pub fn acos(self) -> Self {
    // Based on the Agner Fog "vector class library":
    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
    const_f32_as_f32x4!(P4asinf, 4.2163199048E-2);
    const_f32_as_f32x4!(P3asinf, 2.4181311049E-2);
    const_f32_as_f32x4!(P2asinf, 4.5470025998E-2);
    const_f32_as_f32x4!(P1asinf, 7.4953002686E-2);
    const_f32_as_f32x4!(P0asinf, 1.6666752422E-1);

    let xa = self.abs();
    let big = xa.simd_ge(f32x4::splat(0.5));

    let x1 = f32x4::splat(0.5) * (f32x4::ONE - xa);
    let x2 = xa * xa;
    let x3 = big.blend(x1, x2);

    let xb = x1.sqrt();

    let x4 = big.blend(xb, xa);

    let z = polynomial_4!(x3, P0asinf, P1asinf, P2asinf, P3asinf, P4asinf);
    let z = z.mul_add(x3 * x4, x4);

    let z1 = z + z;

    // acos
    let z3 = self.simd_lt(f32x4::ZERO).blend(f32x4::PI - z1, z1);
    let z4 = f32x4::FRAC_PI_2 - z.flip_signs(self);
    let acos = big.blend(z3, z4);

    acos
  }

  #[inline]
  pub fn atan(self) -> Self {
    // Based on the Agner Fog "vector class library":
    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
    const_f32_as_f32x4!(P3atanf, 8.05374449538E-2);
    const_f32_as_f32x4!(P2atanf, -1.38776856032E-1);
    const_f32_as_f32x4!(P1atanf, 1.99777106478E-1);
    const_f32_as_f32x4!(P0atanf, -3.33329491539E-1);

    let t = self.abs();

    // small:  z = t / 1.0;
    // medium: z = (t-1.0) / (t+1.0);
    // big:    z = -1.0 / t;
    let notsmal = t.simd_ge(Self::SQRT_2 - Self::ONE);
    let notbig = t.simd_le(Self::SQRT_2 + Self::ONE);

    let mut s = notbig.blend(Self::FRAC_PI_4, Self::FRAC_PI_2);
    s = notsmal & s;

    let mut a = notbig & t;
    a = notsmal.blend(a - Self::ONE, a);
    let mut b = notbig & Self::ONE;
    b = notsmal.blend(b + t, b);
    let z = a / b;

    let zz = z * z;

    // Taylor expansion
    let mut re = polynomial_3!(zz, P0atanf, P1atanf, P2atanf, P3atanf);
    re = re.mul_add(zz * z, z) + s;

    // get sign bit
    re = (self.is_sign_negative()).blend(-re, re);

    re
  }

  #[inline]
  pub fn atan2(self, x: Self) -> Self {
    // Based on the Agner Fog "vector class library":
    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
    const_f32_as_f32x4!(P3atanf, 8.05374449538E-2);
    const_f32_as_f32x4!(P2atanf, -1.38776856032E-1);
    const_f32_as_f32x4!(P1atanf, 1.99777106478E-1);
    const_f32_as_f32x4!(P0atanf, -3.33329491539E-1);

    let y = self;

    // move in first octant
    let x1 = x.abs();
    let y1 = y.abs();
    let swapxy = y1.simd_gt(x1);
    // swap x and y if y1 > x1
    let mut x2 = swapxy.blend(y1, x1);
    let mut y2 = swapxy.blend(x1, y1);

    // check for special case: x and y are both +/- INF
    let both_infinite = x.is_inf() & y.is_inf();
    if both_infinite.any() {
      let minus_one = -Self::ONE;
      x2 = both_infinite.blend(x2 & minus_one, x2);
      y2 = both_infinite.blend(y2 & minus_one, y2);
    }

    // x = y = 0 will produce NAN. No problem, fixed below
    let t = y2 / x2;

    // small:  z = t / 1.0;
    // medium: z = (t-1.0) / (t+1.0);
    let notsmal = t.simd_ge(Self::SQRT_2 - Self::ONE);

    let a = notsmal.blend(t - Self::ONE, t);
    let b = notsmal.blend(t + Self::ONE, Self::ONE);
    let s = notsmal & Self::FRAC_PI_4;
    let z = a / b;

    let zz = z * z;

    // Taylor expansion
    let mut re = polynomial_3!(zz, P0atanf, P1atanf, P2atanf, P3atanf);
    re = re.mul_add(zz * z, z) + s;

    // move back in place
    re = swapxy.blend(Self::FRAC_PI_2 - re, re);
    re = ((x | y).simd_eq(Self::ZERO)).blend(Self::ZERO, re);
    re = (x.is_sign_negative()).blend(Self::PI - re, re);

    // get sign bit
    re = (y.is_sign_negative()).blend(-re, re);

    re
  }

  #[inline]
  #[must_use]
  pub fn sin_cos(self) -> (Self, Self) {
    // Based on the Agner Fog "vector class library":
    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h

    const_f32_as_f32x4!(DP1F, 0.78515625_f32 * 2.0);
    const_f32_as_f32x4!(DP2F, 2.4187564849853515625E-4_f32 * 2.0);
    const_f32_as_f32x4!(DP3F, 3.77489497744594108E-8_f32 * 2.0);

    const_f32_as_f32x4!(P0sinf, -1.6666654611E-1);
    const_f32_as_f32x4!(P1sinf, 8.3321608736E-3);
    const_f32_as_f32x4!(P2sinf, -1.9515295891E-4);

    const_f32_as_f32x4!(P0cosf, 4.166664568298827E-2);
    const_f32_as_f32x4!(P1cosf, -1.388731625493765E-3);
    const_f32_as_f32x4!(P2cosf, 2.443315711809948E-5);

    const_f32_as_f32x4!(TWO_OVER_PI, 2.0 / core::f32::consts::PI);

    let xa = self.abs();

    // Find quadrant
    let y = (xa * TWO_OVER_PI).round();
    let q: i32x4 = y.round_int();

    let x = y.mul_neg_add(DP3F, y.mul_neg_add(DP2F, y.mul_neg_add(DP1F, xa)));

    let x2 = x * x;
    let mut s = polynomial_2!(x2, P0sinf, P1sinf, P2sinf) * (x * x2) + x;
    let mut c = polynomial_2!(x2, P0cosf, P1cosf, P2cosf) * (x2 * x2)
      + f32x4::from(0.5).mul_neg_add(x2, f32x4::from(1.0));

    let swap = !(q & i32x4::from(1)).simd_eq(i32x4::from(0));

    let mut overflow: f32x4 = cast(q.simd_gt(i32x4::from(0x2000000)));
    overflow &= xa.is_finite();
    s = overflow.blend(f32x4::from(0.0), s);
    c = overflow.blend(f32x4::from(1.0), c);

    // calc sin
    let mut sin1 = cast::<_, f32x4>(swap).blend(c, s);
    let sign_sin: i32x4 = (q << 30) ^ cast::<_, i32x4>(self);
    sin1 = sin1.flip_signs(cast(sign_sin));

    // calc cos
    let mut cos1 = cast::<_, f32x4>(swap).blend(s, c);
    let sign_cos: i32x4 = ((q + i32x4::from(1)) & i32x4::from(2)) << 30;
    cos1 ^= cast::<_, f32x4>(sign_cos);

    // IEEE 754: sin/cos(±∞) = NaN, sin/cos(NaN) = NaN
    let finite = self.is_finite();
    let nan = Self::splat(f32::NAN);
    let sin_final = finite.blend(sin1, nan);
    let cos_final = finite.blend(cos1, nan);

    (sin_final, cos_final)
  }

  #[inline]
  #[must_use]
  pub fn sin(self) -> Self {
    let (s, _) = self.sin_cos();
    s
  }
  #[inline]
  #[must_use]
  pub fn cos(self) -> Self {
    let (_, c) = self.sin_cos();
    c
  }
  #[inline]
  #[must_use]
  pub fn tan(self) -> Self {
    let (s, c) = self.sin_cos();
    s / c
  }

  /// Calculates hyperbolic sine: `(e^self - e^(-self))/2`.
  #[inline]
  #[must_use]
  pub fn sinh(self) -> Self {
    const_f32_as_f32x4!(P0, 1.0);
    const_f32_as_f32x4!(P1, 1.0 / 6.0);
    const_f32_as_f32x4!(P2, 1.0 / 120.0);
    const_f32_as_f32x4!(P3, 1.0 / 5040.0);
    let a = self.abs();
    // |x| < 0.5: Taylor poly; last truncation term < 1 ULP at x=0.5 for both types
    let small = a.simd_lt(f32x4::from(0.5));
    let t = a * a;
    let poly = a * polynomial_3!(t, P0, P1, P2, P3);
    let exp_based = {
      let e = a.exp();
      (e - Self::ONE / e) * Self::HALF
    };
    let result = small.blend(poly, exp_based);
    result.flip_signs(self)
  }

  /// Calculates hyperbolic cosine: `(e^self + e^(-self))/2`.
  #[inline]
  #[must_use]
  pub fn cosh(self) -> Self {
    const_f32_as_f32x4!(P0, 1.0);
    const_f32_as_f32x4!(P1, 1.0 / 2.0);
    const_f32_as_f32x4!(P2, 1.0 / 24.0);
    const_f32_as_f32x4!(P3, 1.0 / 720.0);
    let a = self.abs();
    // |x| < 0.5: Taylor poly; last truncation term < 1 ULP at x=0.5 for both types
    let small = a.simd_lt(f32x4::from(0.5));
    let t = a * a;
    let poly = polynomial_3!(t, P0, P1, P2, P3);
    let exp_based = {
      let e = a.exp();
      (e + Self::ONE / e) * Self::HALF
    };
    small.blend(poly, exp_based)
  }

  /// Calculates hyperbolic tangent: `sinh(self)/cosh(self)`.
  #[inline]
  #[must_use]
  pub fn tanh(self) -> Self {
    // |x| < 2e-4: tanh(x) ≈ x, error x³/3 < 16·ULP(x)
    // bound: x² < 48·2⁻²³ → x < 2.39e-3; 2e-4 has 10× margin
    // |x| > 9.011: tanh(x) = ±1 to f32 precision (e⁻²ˣ < 2⁻²⁴)
    let a = self.abs();
    let large = a.simd_gt(f32x4::from(9.011));
    if large.all() {
      return Self::ONE.flip_signs(self);
    }
    let small = a.simd_lt(f32x4::from(2e-4));
    let exp_based = {
      let t = (Self::from(-2.0) * a).exp_m1();
      let pos = -t / (t + Self::from(2.0));
      pos.flip_signs(self)
    };
    let result = small.blend(self, exp_based);
    large.blend(Self::ONE.flip_signs(self), result)
  }

  /// Calculates the cube root: `self^(1/3)`.
  #[inline]
  #[must_use]
  pub fn cbrt(self) -> Self {
    let a = self.abs();
    let zero = a.simd_eq(Self::ZERO);
    if zero.all() {
      return self; // preserves -0.0
    }
    let inf = a.is_inf();
    let nan = self.is_nan();

    let tiny = a.simd_lt(Self::from(f32::MIN_POSITIVE));
    let a_work = tiny.blend(a * Self::from(16777216.0), a);

    let e = Self::exponent(a_work) + Self::ONE;
    let d = Self::fraction_2(a_work);

    // C0..C5 from SLEEF's minimax polynomial for 1/cbrt(d) on [0.5, 1.0)
    // Naoki Shibata et al., "SLEEF: A Portable Vectorized Library of C99
    // Mathematical Functions", https://sleef.org / https://github.com/shibatch/sleef
    // Licensed under the Boost Software License 1.0.
    // These are the f32-precision coefficients; our f64 variants use the f64
    // set.
    const_f32_as_f32x4!(C0, 2.2241257);
    const_f32_as_f32x4!(C1, -3.8095417);
    const_f32_as_f32x4!(C2, 5.8982625);
    const_f32_as_f32x4!(C3, -5.532182);
    const_f32_as_f32x4!(C4, 2.8208892);
    const_f32_as_f32x4!(C5, -0.60156447);
    let mut x = polynomial_5!(d, C0, C1, C2, C3, C4, C5);

    let x2 = x * x;
    let x4 = x2 * x2;
    x = x - d.mul_add(x4, -x) * Self::from(1.0 / 3.0);
    // cbrt(d) = d * x² with refinement
    let mut y = (d * x) * x;
    let yx = y * x;
    let t = Self::from(2.0 / 3.0);
    y = y - t * y * (yx - Self::ONE);

    // Scale by 2^(e/3)
    let three = Self::from(3.0);
    let two = Self::from(2.0);
    let neg = e.simd_lt(Self::ZERO);
    let e_adj = neg.blend(e - two, e);
    let k = (e_adj / three).trunc();
    let r = e - three * k;
    const_f32_as_f32x4!(CBRT2, 1.259921);
    const_f32_as_f32x4!(CBRT4, 1.587401);
    y = r.simd_eq(Self::ONE).blend(y * CBRT2, y);
    y = r.simd_eq(two).blend(y * CBRT4, y);
    y *= Self::vm_pow2n(k);
    y = tiny.blend(y / Self::from(256.0_f32), y);

    let result = y.flip_signs(self);
    let result = nan.blend(self, result);
    let result = zero.blend(self, result);
    let result = inf.blend(self, result);
    result
  }

  #[inline]
  #[must_use]
  pub fn to_degrees(self) -> Self {
    const_f32_as_f32x4!(RAD_TO_DEG_RATIO, 180.0_f32 / core::f32::consts::PI);
    self * RAD_TO_DEG_RATIO
  }
  #[inline]
  #[must_use]
  pub fn to_radians(self) -> Self {
    const_f32_as_f32x4!(DEG_TO_RAD_RATIO, core::f32::consts::PI / 180.0_f32);
    self * DEG_TO_RAD_RATIO
  }
  #[inline]
  #[must_use]
  pub fn recip(self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: reciprocal_m128(self.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_div(f32x4_splat(1.0), self.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vdivq_f32(vdupq_n_f32(1.0), self.neon) }}
      } else {
        Self { arr: [
          1.0 / self.arr[0],
          1.0 / self.arr[1],
          1.0 / self.arr[2],
          1.0 / self.arr[3],
        ]}
      }
    }
  }
  #[inline]
  #[must_use]
  pub fn recip_sqrt(self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: reciprocal_sqrt_m128(self.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_div(f32x4_splat(1.0), f32x4_sqrt(self.simd)) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vdivq_f32(vdupq_n_f32(1.0), vsqrtq_f32(self.neon)) }}
      } else if #[cfg(feature="std")] {
        Self { arr: [
          1.0 / self.arr[0].sqrt(),
          1.0 / self.arr[1].sqrt(),
          1.0 / self.arr[2].sqrt(),
          1.0 / self.arr[3].sqrt(),
        ]}
      } else {
        Self { arr: [
          1.0 / software_sqrt(self.arr[0] as f64) as f32,
          1.0 / software_sqrt(self.arr[1] as f64) as f32,
          1.0 / software_sqrt(self.arr[2] as f64) as f32,
          1.0 / software_sqrt(self.arr[3] as f64) as f32,
        ]}
      }
    }
  }
  #[inline]
  #[must_use]
  pub fn sqrt(self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: sqrt_m128(self.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_sqrt(self.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vsqrtq_f32(self.neon) }}
      } else if #[cfg(feature="std")] {
        Self { arr: [
          self.arr[0].sqrt(),
          self.arr[1].sqrt(),
          self.arr[2].sqrt(),
          self.arr[3].sqrt(),
        ]}
      } else {
        Self { arr: [
          software_sqrt(self.arr[0] as f64) as f32,
          software_sqrt(self.arr[1] as f64) as f32,
          software_sqrt(self.arr[2] as f64) as f32,
          software_sqrt(self.arr[3] as f64) as f32,
        ]}
      }
    }
  }

  #[inline]
  #[must_use]
  #[doc(alias("movemask", "move_mask"))]
  pub fn to_bitmask(self) -> u32 {
    pick! {
      if #[cfg(target_feature="sse")] {
        move_mask_m128(self.sse) as u32
      } else if #[cfg(target_feature="simd128")] {
        u32x4_bitmask(self.simd) as u32
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe
        {
          // set all to 1 if top bit is set, else 0
          let masked = vcltq_s32( vreinterpretq_s32_f32(self.neon), vdupq_n_s32(0));

          // select the right bit out of each lane
          let selectbit : uint32x4_t = core::mem::transmute([1u32, 2, 4, 8]);
          let r = vandq_u32(masked, selectbit);

          // horizontally add the 16-bit lanes
          vaddvq_u32(r) as u32
        }
      } else {
        (((self.arr[0].to_bits() as i32) < 0) as u32) << 0 |
        (((self.arr[1].to_bits() as i32) < 0) as u32) << 1 |
        (((self.arr[2].to_bits() as i32) < 0) as u32) << 2 |
        (((self.arr[3].to_bits() as i32) < 0) as u32) << 3
      }
    }
  }
  #[inline]
  #[must_use]
  pub fn any(self) -> bool {
    pick! {
      if #[cfg(target_feature="simd128")] {
        v128_any_true(self.simd)
      } else {
        self.to_bitmask() != 0
      }
    }
  }
  #[inline]
  #[must_use]
  pub fn all(self) -> bool {
    pick! {
      if #[cfg(target_feature="simd128")] {
        u32x4_all_true(self.simd)
      } else {
        // four lanes
        self.to_bitmask() == 0b1111
      }
    }
  }
  #[inline]
  #[must_use]
  pub fn none(self) -> bool {
    !self.any()
  }

  #[inline]
  fn vm_pow2n(self) -> Self {
    const_f32_as_f32x4!(pow2_23, 8388608.0);
    const_f32_as_f32x4!(bias, 127.0);
    let a = self + (bias + pow2_23);
    let c = cast::<_, i32x4>(a) << 23;
    let std_result = cast::<_, f32x4>(c);

    let min_exp = f32x4::from(-126.0);
    let is_sub = self.simd_lt(min_exp);
    if is_sub.any() {
      let valid = self.simd_ge(f32x4::from(-149.0));
      let shift_f = self + f32x4::from(149.0);
      let mut shift_i = shift_f.trunc_int();
      shift_i = cast::<_, i32x4>(valid).blend(shift_i, i32x4::ZERO);
      let mantissa = i32x4::ONE << shift_i;
      let sub_result = cast::<_, f32x4>(mantissa);
      let sub_result = valid.blend(sub_result, f32x4::ZERO);
      is_sub.blend(sub_result, std_result)
    } else {
      std_result
    }
  }

  /// Calculate the exponent of a packed `f32x4`
  #[inline]
  #[must_use]
  pub fn exp(self) -> Self {
    const_f32_as_f32x4!(P0, 1.0 / 2.0);
    const_f32_as_f32x4!(P1, 1.0 / 6.0);
    const_f32_as_f32x4!(P2, 1.0 / 24.0);
    const_f32_as_f32x4!(P3, 1.0 / 120.0);
    const_f32_as_f32x4!(P4, 1.0 / 720.0);
    const_f32_as_f32x4!(P5, 1.0 / 5040.0);
    // LN2D_HI/LO: double-double decomposition of ln(2) for exp range reduction,
    // following the approach from fdlibm's e_exp.c (Sun Microsystems,
    // https://www.netlib.org/fdlibm/). The f32 split uses f32-precision constants
    // (0.693359375, -2.12194440e-4) summing to ln(2) with single-precision
    // accuracy; the f64 variants use a full f64 double-double
    // decomposition.
    const_f32_as_f32x4!(LN2D_HI, 0.693359375);
    const_f32_as_f32x4!(LN2D_LO, -2.12194440e-4);
    // max_x = ln(f32::MAX) ≈ 88.7229, max_r = 127 (IEEE max normal exponent)
    // min_x = -149.5 ln(2) ≈ -103.63: min r for vm_pow2n subnormal
    let max_x = f32x4::from(88.723);
    let min_x = f32x4::from(-103.63);
    // x < min_x: e^x underflows to 0 -- skip the entire pipeline
    let finite = self.is_finite();
    let neg_underflow = self.simd_lt(min_x) & finite;
    if neg_underflow.all() {
      return Self::ZERO;
    }
    let max_r = f32x4::from(127.0);
    let r = (self * Self::LOG2_E).round();
    let big = r.simd_gt(max_r);
    let r_safe = big.blend(max_r, r);
    let excess = r - max_r;
    let excess = big.blend(excess, Self::ZERO);
    let scale = Self::vm_pow2n(excess);
    let x = r.mul_neg_add(LN2D_HI, self);
    let x = r.mul_neg_add(LN2D_LO, x);
    let z = polynomial_5!(x, P0, P1, P2, P3, P4, P5);
    let x2 = x * x;
    let z = z.mul_add(x2, x);
    let n2 = Self::vm_pow2n(r_safe);
    let z = (z + Self::ONE) * scale * n2;
    let nan_mask = self.is_nan();
    let mut result = nan_mask.blend(Self::nan_pow(), z);
    let pos_overflow = self.simd_gt(max_x) & finite;
    result = pos_overflow.blend(Self::infinity(), result);
    result = neg_underflow.blend(Self::ZERO, result);
    let pos_inf = !finite & !self.is_sign_negative() & !nan_mask;
    result = pos_inf.blend(Self::infinity(), result);
    let neg_inf = !finite & self.is_sign_negative() & !nan_mask;
    result = neg_inf.blend(Self::ZERO, result);
    result
  }

  /// Calculate `e^self - 1` for each lane.
  /// Accurate even for very small values.
  #[inline]
  #[must_use]
  pub fn exp_m1(self) -> Self {
    // x < -17.329: e^x < 2⁻²⁵, exp_m1(x) = -1.0 exactly (mantissa exhaustion)
    // IEEE simd_lt returns false for NaN, so NaN lanes can't reach here.
    // -inf is < -17.329, and exp_m1(-inf) = -1.0, also correct.
    if self.simd_lt(f32x4::from(-17.329)).all() {
      return f32x4::from(-1.0);
    }
    const_f32_as_f32x4!(P0, 1.0 / 2.0);
    const_f32_as_f32x4!(P1, 1.0 / 6.0);
    const_f32_as_f32x4!(P2, 1.0 / 24.0);
    const_f32_as_f32x4!(P3, 1.0 / 120.0);
    const_f32_as_f32x4!(P4, 1.0 / 720.0);
    const_f32_as_f32x4!(P5, 1.0 / 5040.0);
    // LN2D_HI/LO: double-double decomposition of ln(2) for exp range reduction,
    // following the approach from fdlibm's e_exp.c (Sun Microsystems,
    // https://www.netlib.org/fdlibm/). The f32 split uses f32-precision constants
    // (0.693359375, -2.12194440e-4) summing to ln(2) with single-precision
    // accuracy; the f64 variants use a full f64 double-double
    // decomposition.
    const_f32_as_f32x4!(LN2D_HI, 0.693359375);
    const_f32_as_f32x4!(LN2D_LO, -2.12194440e-4);
    // max_x = ln(f32::MAX) ≈ 88.7229, max_r = 127 (IEEE max normal exponent)
    // min_x = -149.5 ln(2) ≈ -103.63: min r for vm_pow2n subnormal
    let max_x = f32x4::from(88.723);
    let min_x = f32x4::from(-103.63);
    let max_r = f32x4::from(127.0);
    let r = (self * Self::LOG2_E).round();
    let big = r.simd_gt(max_r);
    let r_safe = big.blend(max_r, r);
    let excess = r - max_r;
    let excess = big.blend(excess, Self::ZERO);
    let scale = Self::vm_pow2n(excess);
    let x = r.mul_neg_add(LN2D_HI, self);
    let x = r.mul_neg_add(LN2D_LO, x);
    let z = polynomial_5!(x, P0, P1, P2, P3, P4, P5);
    let x2 = x * x;
    let z = z.mul_add(x2, x);
    let n2 = Self::vm_pow2n(r_safe);
    let exp_val = (z + Self::ONE) * scale * n2;
    let r_is_zero = r.simd_eq(Self::ZERO);
    let z = r_is_zero.blend(z, exp_val - Self::ONE);
    let nan_mask = self.is_nan();
    let finite = self.is_finite();
    let mut result = nan_mask.blend(Self::nan_pow(), z);
    let pos_overflow = self.simd_gt(max_x) & finite;
    result = pos_overflow.blend(Self::infinity(), result);
    let neg_underflow = self.simd_lt(min_x) & finite;
    result = neg_underflow.blend(-Self::ONE, result);
    let pos_inf = !finite & !self.is_sign_negative() & !nan_mask;
    result = pos_inf.blend(Self::infinity(), result);
    let neg_inf = !finite & self.is_sign_negative() & !nan_mask;
    result = neg_inf.blend(-Self::ONE, result);
    let is_zero = self.simd_eq(Self::ZERO);
    result = is_zero.blend(self, result);
    result
  }

  /// Returns `2^self`.
  #[inline]
  #[must_use]
  pub fn exp2(self) -> Self {
    const_f32_as_f32x4!(P2, 1.0 / 2.0);
    const_f32_as_f32x4!(P3, 1.0 / 6.0);
    const_f32_as_f32x4!(P4, 1.0 / 24.0);
    const_f32_as_f32x4!(P5, 1.0 / 120.0);
    const_f32_as_f32x4!(P6, 1.0 / 720.0);
    const_f32_as_f32x4!(P7, 1.0 / 5040.0);

    // max_x = log2(f32::MAX) ≈ 127.99999
    // min_x = log2(f32::MIN_POSITIVE) - 23 ≈ -126 - 23 = -149
    let max_x = f32x4::from(127.99999);
    let min_x = f32x4::from(-149.5);
    let finite = self.is_finite();
    let neg_underflow = self.simd_lt(min_x) & finite;
    if neg_underflow.all() {
      return Self::ZERO;
    }

    let round = self.round();
    let max_r = f32x4::from(127.0);
    let big = round.simd_gt(max_r);
    let r_safe = big.blend(max_r, round);
    let excess = round - max_r;
    let excess = big.blend(excess, Self::ZERO);
    let scale = Self::vm_pow2n(excess);

    let fract = (self - round) * Self::LN_2;
    let fract_partial_exp2 = polynomial_5!(fract, P2, P3, P4, P5, P6, P7);
    let fract2 = fract * fract;
    let fract_exp2 = fract_partial_exp2.mul_add(fract2, fract) + Self::ONE;

    let n2 = Self::vm_pow2n(r_safe);
    let result = fract_exp2 * scale * n2;

    let nan_mask = self.is_nan();
    let mut result = nan_mask.blend(Self::nan_pow(), result);
    let pos_overflow = self.simd_gt(max_x) & finite;
    result = pos_overflow.blend(Self::infinity(), result);
    result = neg_underflow.blend(Self::ZERO, result);
    let pos_inf = !finite & !self.is_sign_negative() & !nan_mask;
    result = pos_inf.blend(Self::infinity(), result);
    let neg_inf = !finite & self.is_sign_negative() & !nan_mask;
    result = neg_inf.blend(Self::ZERO, result);
    result
  }

  #[inline]
  fn exponent(self) -> f32x4 {
    const_f32_as_f32x4!(pow2_23, 8388608.0);
    const_f32_as_f32x4!(bias, 127.0);
    let a = cast::<_, u32x4>(self);
    let b = a >> 23;
    let c = b | cast::<_, u32x4>(pow2_23);
    let d = cast::<_, f32x4>(c);
    let e = d - (pow2_23 + bias);
    e
  }

  #[inline]
  fn fraction_2(self) -> Self {
    let t1 = cast::<_, u32x4>(self);
    let t2 = cast::<_, u32x4>(
      (t1 & u32x4::from(0x007FFFFF)) | u32x4::from(0x3F000000),
    );
    cast::<_, f32x4>(t2)
  }
  #[inline]
  fn is_zero_or_subnormal(self) -> Self {
    let t = cast::<_, i32x4>(self);
    let t = t & i32x4::splat(0x7F800000);
    let mask = t.simd_eq(i32x4::splat(0));
    cast::<_, f32x4>(mask)
  }
  #[inline]
  fn infinity() -> Self {
    cast::<_, f32x4>(i32x4::splat(0x7F800000))
  }
  #[inline]
  fn nan_log() -> Self {
    cast::<_, f32x4>(i32x4::splat(0x7FC00000 | 0x101 & 0x003FFFFF))
  }
  #[inline]
  fn nan_pow() -> Self {
    cast::<_, f32x4>(i32x4::splat(0x7FC00000 | 0x101 & 0x003FFFFF))
  }

  /// Returns true for each element if it has a positive sign, including `+0.0`,
  /// `NaN`s with positive sign bit and positive infinity.
  #[inline]
  #[must_use]
  pub fn is_sign_positive(self) -> Self {
    const SIGN_MASK: u32x4 = u32x4::splat((-0.0_f32).to_bits());

    let bits = cast::<f32x4, u32x4>(self);
    let sign = bits & SIGN_MASK;
    let result = sign.simd_eq(u32x4::ZERO);
    cast::<u32x4, f32x4>(result)
  }

  /// Returns true for each element if it has a negative sign, including `-0.0`,
  /// `NaN`s with negative sign bit and negative infinity.
  #[inline]
  #[must_use]
  pub fn is_sign_negative(self) -> Self {
    const SIGN_MASK: u32x4 = u32x4::splat((-0.0_f32).to_bits());

    let bits = cast::<f32x4, u32x4>(self);
    let sign = bits & SIGN_MASK;
    let result = sign.simd_eq(SIGN_MASK);
    cast::<u32x4, f32x4>(result)
  }

  /// horizontal add of all the elements of the vector
  #[inline]
  #[must_use]
  pub fn reduce_add(self) -> f32 {
    let arr: [f32; 4] = cast(self);
    arr.iter().sum()
  }

  /// horizontal multiplication of all the elements of the vector
  #[inline]
  #[must_use]
  pub fn reduce_mul(self) -> f32 {
    let arr: [f32; 4] = cast(self);
    arr.iter().product()
  }

  /// Natural log (ln(x))
  #[inline]
  #[must_use]
  pub fn ln(self) -> Self {
    const_f32_as_f32x4!(HALF, 0.5);
    const_f32_as_f32x4!(P0, 3.3333331174E-1);
    const_f32_as_f32x4!(P1, -2.4999993993E-1);
    const_f32_as_f32x4!(P2, 2.0000714765E-1);
    const_f32_as_f32x4!(P3, -1.6668057665E-1);
    const_f32_as_f32x4!(P4, 1.4249322787E-1);
    const_f32_as_f32x4!(P5, -1.2420140846E-1);
    const_f32_as_f32x4!(P6, 1.1676998740E-1);
    const_f32_as_f32x4!(P7, -1.1514610310E-1);
    const_f32_as_f32x4!(P8, 7.0376836292E-2);
    const_f32_as_f32x4!(LN2F_HI, 0.693359375);
    const_f32_as_f32x4!(LN2F_LO, -2.12194440e-4);
    const_f32_as_f32x4!(VM_SMALLEST_NORMAL, 1.17549435E-38);

    let x1 = self;
    let x = Self::fraction_2(x1);
    let e = Self::exponent(x1);
    let mask = x.simd_gt(Self::SQRT_2 * HALF);
    let x = (!mask).blend(x + x, x);
    let fe = mask.blend(e + Self::ONE, e);
    let x = x - Self::ONE;
    let res = polynomial_8!(x, P0, P1, P2, P3, P4, P5, P6, P7, P8);
    let x2 = x * x;
    let res = x2 * x * res;
    let res = fe.mul_add(LN2F_LO, res);
    let res = res + x2.mul_neg_add(HALF, x);
    let res = fe.mul_add(LN2F_HI, res);
    let overflow = !self.is_finite();
    let underflow = x1.simd_lt(VM_SMALLEST_NORMAL);
    let mask = overflow | underflow;
    if !mask.any() {
      res
    } else {
      let is_zero = self.is_zero_or_subnormal();
      let res = underflow.blend(Self::nan_log(), res);
      // Note: is_zero_or_subnormal() lumps subnormals (exponent==0) with zero.
      // Both get -Inf here. True subnormal inputs (~1.4e-45..1.175e-38) should
      // produce a finite negative result, but are vanishingly rare in
      // practice.
      let res = is_zero.blend(-Self::infinity(), res);
      let res = overflow.blend(self, res);
      // This must come *after* overflow.blend to overwrite ln(-∞) = -∞ to NaN
      let res = (!self.is_finite() & self.is_sign_negative())
        .blend(Self::nan_log(), res);
      res
    }
  }

  /// Calculate `ln(1 + self)` for each lane.
  /// Accurate even for very small values.
  #[inline]
  #[must_use]
  pub fn ln_1p(self) -> Self {
    // Based on the identity ln(1+x) = x·ln(1+x)/((1+x)-1), i.e. x·ln(u)/(u-1)
    // where u = 1+x. From MUSL libc (Rich Felker et al., https://musl.libc.org) src/math/log1pf.c
    // and fdlibm (Sun Microsystems, https://www.netlib.org/fdlibm/) s_log1p.c.
    // When 1+x rounds to 1 exactly (subnormal x), return x directly.
    // When 1+x overflows (+inf), return ln(u) without correction.
    // Mathematically exact: compensates for the rounding loss in 1+x without
    // needing a series threshold.
    let u = self + Self::ONE;
    let eq = u.simd_eq(Self::ONE);
    let ln_u = Self::ln(u);
    let correction = self * (ln_u / (u - Self::ONE));
    let result = eq.blend(self, correction);
    let over = u.is_inf();
    over.blend(ln_u, result)
  }

  #[inline]
  #[must_use]
  pub fn log2(self) -> Self {
    Self::ln(self) * Self::LOG2_E
  }
  #[inline]
  #[must_use]
  pub fn log10(self) -> Self {
    Self::ln(self) * Self::LOG10_E
  }

  #[inline]
  #[must_use]
  pub fn pow_f32x4(self, y: f32x4) -> Self {
    const_f32_as_f32x4!(ln2f_hi, 0.693359375);
    const_f32_as_f32x4!(ln2f_lo, -2.12194440e-4);
    const_f32_as_f32x4!(P0logf, 3.3333331174E-1);
    const_f32_as_f32x4!(P1logf, -2.4999993993E-1);
    const_f32_as_f32x4!(P2logf, 2.0000714765E-1);
    const_f32_as_f32x4!(P3logf, -1.6668057665E-1);
    const_f32_as_f32x4!(P4logf, 1.4249322787E-1);
    const_f32_as_f32x4!(P5logf, -1.2420140846E-1);
    const_f32_as_f32x4!(P6logf, 1.1676998740E-1);
    const_f32_as_f32x4!(P7logf, -1.1514610310E-1);
    const_f32_as_f32x4!(P8logf, 7.0376836292E-2);

    const_f32_as_f32x4!(p2expf, 1.0 / 2.0); // coefficients for Taylor expansion of exp
    const_f32_as_f32x4!(p3expf, 1.0 / 6.0);
    const_f32_as_f32x4!(p4expf, 1.0 / 24.0);
    const_f32_as_f32x4!(p5expf, 1.0 / 120.0);
    const_f32_as_f32x4!(p6expf, 1.0 / 720.0);
    const_f32_as_f32x4!(p7expf, 1.0 / 5040.0);

    let x1 = self.abs();
    let x = x1.fraction_2();

    let mask = x.simd_gt(f32x4::SQRT_2 * f32x4::HALF);
    let x = (!mask).blend(x + x, x);

    let x = x - f32x4::ONE;
    let x2 = x * x;
    let lg1 = polynomial_8!(
      x, P0logf, P1logf, P2logf, P3logf, P4logf, P5logf, P6logf, P7logf, P8logf
    );
    let lg1 = lg1 * x2 * x;

    let ef = x1.exponent();
    let ef = mask.blend(ef + f32x4::ONE, ef);

    let e1 = (ef * y).round();
    let yr = ef.mul_sub(y, e1);

    let lg = f32x4::HALF.mul_neg_add(x2, x) + lg1;
    let x2_err = (f32x4::HALF * x).mul_sub(x, f32x4::HALF * x2);
    let lg_err = f32x4::HALF.mul_add(x2, lg - x) - lg1;

    let e2 = (lg * y * f32x4::LOG2_E).round();
    let v = lg.mul_sub(y, e2 * ln2f_hi);
    let v = e2.mul_neg_add(ln2f_lo, v);
    let v = v - (lg_err + x2_err).mul_sub(y, yr * f32x4::LN_2);

    let x = v;
    let e3 = (x * f32x4::LOG2_E).round();
    let x = e3.mul_neg_add(f32x4::LN_2, x);
    let x2 = x * x;
    let z = x2.mul_add(
      polynomial_5!(x, p2expf, p3expf, p4expf, p5expf, p6expf, p7expf),
      x + f32x4::ONE,
    );

    let ee = e1 + e2 + e3;
    let ei = cast::<_, i32x4>(ee.round_int());
    let ej = cast::<_, i32x4>(ei + (cast::<_, i32x4>(z) >> 23));

    let overflow = cast::<_, f32x4>(ej.simd_gt(i32x4::splat(0x0FF)))
      | (ee.simd_gt(f32x4::splat(300.0)));
    let underflow = cast::<_, f32x4>(ej.simd_lt(i32x4::splat(0x000)))
      | (ee.simd_lt(f32x4::splat(-300.0)));

    // Add exponent by integer addition
    let z = cast::<_, f32x4>(cast::<_, i32x4>(z) + (ei << 23));

    // Check for overflow/underflow
    let z = if (overflow | underflow).any() {
      let z = underflow.blend(f32x4::ZERO, z);
      overflow.blend(Self::infinity(), z)
    } else {
      z
    };

    // Check for self == 0
    let x_zero = self.is_zero_or_subnormal();
    let z = x_zero.blend(
      y.simd_lt(f32x4::ZERO).blend(
        Self::infinity(),
        y.simd_eq(f32x4::ZERO).blend(f32x4::ONE, f32x4::ZERO),
      ),
      z,
    );

    let x_sign = self.is_sign_negative();
    let z = if x_sign.any() {
      // Y into an integer
      let yi = y.simd_eq(y.round());
      // Is y odd?
      let y_odd = cast::<_, i32x4>(y.round_int() << 31).round_float();

      let z1 =
        yi.blend(z | y_odd, self.simd_eq(Self::ZERO).blend(z, Self::nan_pow()));
      x_sign.blend(z1, z)
    } else {
      z
    };

    let x_finite = self.is_finite();
    let y_finite = y.is_finite();
    let e_finite = ee.is_finite();
    if (x_finite & y_finite & (e_finite | x_zero)).all() {
      return z;
    }

    (self.is_nan() | y.is_nan()).blend(self + y, z)
  }

  #[inline]
  pub fn powf(self, y: f32) -> Self {
    Self::pow_f32x4(self, f32x4::splat(y))
  }

  #[must_use]
  #[inline]
  pub fn unpack_lo(self, b: Self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: unpack_low_m128(self.sse, b.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self {
          simd: u32x4_shuffle::<0, 4, 1, 5>(self.simd, b.simd)
        }
      } else if #[cfg(all(target_feature="neon", target_arch="aarch64"))]{
        unsafe {Self { neon: vzip1q_f32(self.neon, b.neon) }}
      } else {
        Self { arr: [
          self.arr[0],
          b.arr[0],
          self.arr[1],
          b.arr[1],
        ]}
      }
    }
  }

  #[must_use]
  #[inline]
  pub fn unpack_hi(self, b: Self) -> Self {
    pick! {
      if #[cfg(target_feature="sse")] {
        Self { sse: unpack_high_m128(self.sse, b.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self {
          simd: u32x4_shuffle::<2, 6, 3, 7>(self.simd, b.simd)
        }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
        unsafe {Self { neon: vzip2q_f32(self.neon, b.neon) }}
      } else {
        Self { arr: [
          self.arr[2],
          b.arr[2],
          self.arr[3],
          b.arr[3],
        ]}
      }
    }
  }

  /// Transpose matrix of 4x4 `f32` matrix. Currently only accelerated on SSE.
  #[must_use]
  #[inline]
  pub fn transpose(data: [f32x4; 4]) -> [f32x4; 4] {
    pick! {
      if #[cfg(target_feature="sse")] {
        let mut e0 = data[0];
        let mut e1 = data[1];
        let mut e2 = data[2];
        let mut e3 = data[3];

        transpose_four_m128(&mut e0.sse, &mut e1.sse, &mut e2.sse, &mut e3.sse);

        [e0, e1, e2, e3]
      } else if #[cfg(any(all(target_feature="neon",target_arch="aarch64"), target_feature="simd128"))] {
        let a = data[0].unpack_lo(data[2]);
        let b = data[1].unpack_lo(data[3]);
        let c = data[0].unpack_hi(data[2]);
        let d = data[1].unpack_hi(data[3]);

        [
          a.unpack_lo(b),
          a.unpack_hi(b),
          c.unpack_lo(d),
          c.unpack_hi(d),
        ]
      } else {
        #[inline(always)]
        fn transpose_column(data: &[f32x4; 4], index: usize) -> f32x4 {
          f32x4::new([
            data[0].as_array()[index],
            data[1].as_array()[index],
            data[2].as_array()[index],
            data[3].as_array()[index],
          ])
        }

        [
          transpose_column(&data, 0),
          transpose_column(&data, 1),
          transpose_column(&data, 2),
          transpose_column(&data, 3),
        ]
      }
    }
  }

  #[inline]
  pub fn to_array(self) -> [f32; 4] {
    cast(self)
  }

  #[inline]
  pub fn as_array(&self) -> &[f32; 4] {
    cast_ref(self)
  }

  #[inline]
  pub fn as_mut_array(&mut self) -> &mut [f32; 4] {
    cast_mut(self)
  }

  #[inline]
  pub fn from_i32x4(v: i32x4) -> Self {
    pick! {
      if #[cfg(target_feature="sse2")] {
        Self { sse: convert_to_m128_from_i32_m128i(v.sse) }
      } else if #[cfg(target_feature="simd128")] {
        Self { simd: f32x4_convert_i32x4(v.simd) }
      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
        Self { neon: unsafe { vcvtq_f32_s32(v.neon) }}
      } else {
        Self { arr: [
            v.as_array()[0] as f32,
            v.as_array()[1] as f32,
            v.as_array()[2] as f32,
            v.as_array()[3] as f32,
          ] }
      }
    }
  }

  /// Returns true for each element if its sign bit is set.
  ///
  /// If the sign bit is set, the result has all bits set, not just the sign
  /// bit. This has been renamed to [`is_sign_negative`].
  ///
  /// [`is_sign_negative`]: Self::is_sign_negative
  #[inline]
  #[must_use]
  #[deprecated(since = "1.4.0", note = "renamed to `is_sign_negative`")]
  pub fn sign_bit(self) -> Self {
    self.is_sign_negative()
  }
}