wide/
f32x4_.rs

1use super::*;
2
3pick! {
4  if #[cfg(target_feature="sse")] {
5    /// A SIMD vector with four elements of type [`f32`].
6    ///
7    /// See the [crate level documentation] for more information about SIMD
8    /// vectors.
9    ///
10    /// [crate level documentation]: crate
11    #[derive(Default, Clone, Copy, PartialEq)]
12    #[repr(C, align(16))]
13    pub struct f32x4 { pub(crate) sse: m128 }
14  } else if #[cfg(target_feature="simd128")] {
15    use core::arch::wasm32::*;
16
17    /// A SIMD vector with four elements of type [`f32`].
18    ///
19    /// See the [crate level documentation] for more information about SIMD
20    /// vectors.
21    ///
22    /// [crate level documentation]: crate
23    #[derive(Clone, Copy)]
24    #[repr(transparent)]
25    pub struct f32x4 { pub(crate) simd: v128 }
26
27    impl Default for f32x4 {
28      fn default() -> Self {
29        Self::splat(0.0)
30      }
31    }
32
33    impl PartialEq for f32x4 {
34      fn eq(&self, other: &Self) -> bool {
35        u32x4_all_true(f32x4_eq(self.simd, other.simd))
36      }
37    }
38  } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
39    use core::arch::aarch64::*;
40
41    /// A SIMD vector with four elements of type [`f32`].
42    ///
43    /// See the [crate level documentation] for more information about SIMD
44    /// vectors.
45    ///
46    /// [crate level documentation]: crate
47    #[repr(C)]
48    #[derive(Copy, Clone)]
49    pub struct f32x4 { pub(crate) neon : float32x4_t }
50
51    impl Default for f32x4 {
52      #[inline]
53      fn default() -> Self {
54        unsafe { Self { neon: vdupq_n_f32(0.0)} }
55      }
56    }
57
58    impl PartialEq for f32x4 {
59      #[inline]
60      fn eq(&self, other: &Self) -> bool {
61        unsafe { vminvq_u32(vceqq_f32(self.neon, other.neon))==u32::MAX }
62      }
63
64    }
65    } else {
66    /// A SIMD vector with four elements of type [`f32`].
67    ///
68    /// See the [crate level documentation] for more information about SIMD
69    /// vectors.
70    ///
71    /// [crate level documentation]: crate
72    #[derive(Default, Clone, Copy, PartialEq)]
73    #[repr(C, align(16))]
74    pub struct f32x4 { pub(crate) arr: [f32;4] }
75  }
76}
77
78macro_rules! const_f32_as_f32x4 {
79  ($i:ident, $f:expr) => {
80    #[allow(non_upper_case_globals)]
81    pub const $i: f32x4 = f32x4::new([$f; 4]);
82  };
83}
84
85impl_simd_float! {
86  unsafe {
87    T = f32,
88    N = 4,
89    Simd = f32x4,
90    IntSimd = i32x4,
91    UintT = u32,
92    UintSimd = u32x4,
93    optional_type_x86_inner { X86Inner = __m128 },
94    optional_type_arm_inner { ArmInner = float32x4_t },
95    optional_type_wasm_inner { WasmInner = v128 },
96  }
97  old_powf_simd_fn_name = pow_f32x4,
98
99  #[inline]
100  fn neg(self) -> Self::Output {
101    pick! {
102      if #[cfg(target_feature="sse")] {
103        Self { sse: bitxor_m128(self.sse, Self::splat(-0.0).sse) }
104      } else if #[cfg(target_feature="simd128")] {
105        Self { simd: f32x4_neg(self.simd) }
106      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
107        unsafe {Self { neon: vnegq_f32(self.neon) }}
108      } else {
109        Self { arr: [
110          -self.arr[0],
111          -self.arr[1],
112          -self.arr[2],
113          -self.arr[3],
114        ]}
115      }
116    }
117  }
118
119  #[inline]
120  fn not(self) -> Self::Output {
121    self ^ cast::<u128, f32x4>(u128::MAX)
122  }
123
124  #[inline]
125  fn add(self, rhs: Self) -> Self::Output {
126    pick! {
127      if #[cfg(target_feature="sse")] {
128        Self { sse: add_m128(self.sse, rhs.sse) }
129      } else if #[cfg(target_feature="simd128")] {
130        Self { simd: f32x4_add(self.simd, rhs.simd) }
131      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
132        unsafe { Self { neon: vaddq_f32(self.neon, rhs.neon) } }
133      } else {
134        Self { arr: [
135          self.arr[0] + rhs.arr[0],
136          self.arr[1] + rhs.arr[1],
137          self.arr[2] + rhs.arr[2],
138          self.arr[3] + rhs.arr[3],
139        ]}
140      }
141    }
142  }
143
144  #[inline]
145  fn sub(self, rhs: Self) -> Self::Output {
146    pick! {
147      if #[cfg(target_feature="sse")] {
148        Self { sse: sub_m128(self.sse, rhs.sse) }
149      } else if #[cfg(target_feature="simd128")] {
150        Self { simd: f32x4_sub(self.simd, rhs.simd) }
151      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
152        unsafe {Self { neon: vsubq_f32(self.neon, rhs.neon) }}
153      } else {
154        Self { arr: [
155          self.arr[0] - rhs.arr[0],
156          self.arr[1] - rhs.arr[1],
157          self.arr[2] - rhs.arr[2],
158          self.arr[3] - rhs.arr[3],
159        ]}
160      }
161    }
162  }
163
164  #[inline]
165  fn mul(self, rhs: Self) -> Self::Output {
166    pick! {
167      if #[cfg(target_feature="sse")] {
168        Self { sse: mul_m128(self.sse, rhs.sse) }
169      } else if #[cfg(target_feature="simd128")] {
170        Self { simd: f32x4_mul(self.simd, rhs.simd) }
171      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
172        unsafe {Self { neon: vmulq_f32(self.neon, rhs.neon) }}
173      } else {
174        Self { arr: [
175          self.arr[0] * rhs.arr[0],
176          self.arr[1] * rhs.arr[1],
177          self.arr[2] * rhs.arr[2],
178          self.arr[3] * rhs.arr[3],
179        ]}
180      }
181    }
182  }
183
184  #[inline]
185  fn div(self, rhs: Self) -> Self::Output {
186    pick! {
187      if #[cfg(target_feature="sse")] {
188        Self { sse: div_m128(self.sse, rhs.sse) }
189      } else if #[cfg(target_feature="simd128")] {
190        Self { simd: f32x4_div(self.simd, rhs.simd) }
191      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
192        unsafe {Self { neon: vdivq_f32(self.neon, rhs.neon) }}
193      } else {
194        Self { arr: [
195          self.arr[0] / rhs.arr[0],
196          self.arr[1] / rhs.arr[1],
197          self.arr[2] / rhs.arr[2],
198          self.arr[3] / rhs.arr[3],
199        ]}
200      }
201    }
202  }
203
204  #[inline]
205  fn rem(self, rhs: Self) -> Self::Output {
206    Self::new([
207      self.to_array()[0] % rhs.to_array()[0],
208      self.to_array()[1] % rhs.to_array()[1],
209      self.to_array()[2] % rhs.to_array()[2],
210      self.to_array()[3] % rhs.to_array()[3],
211    ])
212  }
213
214  #[inline]
215  fn bitand(self, rhs: Self) -> Self::Output {
216    pick! {
217      if #[cfg(target_feature="sse")] {
218        Self { sse: bitand_m128(self.sse, rhs.sse) }
219      } else if #[cfg(target_feature="simd128")] {
220        Self { simd: v128_and(self.simd, rhs.simd) }
221      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
222        unsafe {Self { neon: vreinterpretq_f32_u32(vandq_u32(vreinterpretq_u32_f32(self.neon), vreinterpretq_u32_f32(rhs.neon))) }}
223      } else {
224        Self { arr: [
225          f32::from_bits(self.arr[0].to_bits() & rhs.arr[0].to_bits()),
226          f32::from_bits(self.arr[1].to_bits() & rhs.arr[1].to_bits()),
227          f32::from_bits(self.arr[2].to_bits() & rhs.arr[2].to_bits()),
228          f32::from_bits(self.arr[3].to_bits() & rhs.arr[3].to_bits()),
229        ]}
230      }
231    }
232  }
233
234  #[inline]
235  fn bitor(self, rhs: Self) -> Self::Output {
236    pick! {
237      if #[cfg(target_feature="sse")] {
238        Self { sse: bitor_m128(self.sse, rhs.sse) }
239      } else if #[cfg(target_feature="simd128")] {
240        Self { simd: v128_or(self.simd, rhs.simd) }
241      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
242        unsafe {Self { neon: vreinterpretq_f32_u32(vorrq_u32(vreinterpretq_u32_f32(self.neon), vreinterpretq_u32_f32(rhs.neon))) }}
243      } else {
244        Self { arr: [
245          f32::from_bits(self.arr[0].to_bits() | rhs.arr[0].to_bits()),
246          f32::from_bits(self.arr[1].to_bits() | rhs.arr[1].to_bits()),
247          f32::from_bits(self.arr[2].to_bits() | rhs.arr[2].to_bits()),
248          f32::from_bits(self.arr[3].to_bits() | rhs.arr[3].to_bits()),
249        ]}
250      }
251    }
252  }
253
254  #[inline]
255  fn bitxor(self, rhs: Self) -> Self::Output {
256    pick! {
257      if #[cfg(target_feature="sse")] {
258        Self { sse: bitxor_m128(self.sse, rhs.sse) }
259      } else if #[cfg(target_feature="simd128")] {
260        Self { simd: v128_xor(self.simd, rhs.simd) }
261      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
262        unsafe {Self { neon: vreinterpretq_f32_u32(veorq_u32(vreinterpretq_u32_f32(self.neon), vreinterpretq_u32_f32(rhs.neon))) }}
263      } else {
264        Self { arr: [
265          f32::from_bits(self.arr[0].to_bits() ^ rhs.arr[0].to_bits()),
266          f32::from_bits(self.arr[1].to_bits() ^ rhs.arr[1].to_bits()),
267          f32::from_bits(self.arr[2].to_bits() ^ rhs.arr[2].to_bits()),
268          f32::from_bits(self.arr[3].to_bits() ^ rhs.arr[3].to_bits()),
269        ]}
270      }
271    }
272  }
273
274  #[inline]
275  fn simd_eq(self, rhs: Self) -> Self::Output {
276    pick! {
277      if #[cfg(target_feature="sse")] {
278        Self { sse: cmp_eq_mask_m128(self.sse, rhs.sse) }
279      } else if #[cfg(target_feature="simd128")] {
280        Self { simd: f32x4_eq(self.simd, rhs.simd) }
281      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
282        unsafe {Self { neon: vreinterpretq_f32_u32(vceqq_f32(self.neon, rhs.neon)) }}
283      } else {
284        Self { arr: [
285          if self.arr[0] == rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
286          if self.arr[1] == rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
287          if self.arr[2] == rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
288          if self.arr[3] == rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
289        ]}
290      }
291    }
292  }
293
294  #[inline]
295  fn simd_ne(self, rhs: Self) -> Self::Output {
296    pick! {
297      if #[cfg(target_feature="sse")] {
298        Self { sse: cmp_neq_mask_m128(self.sse, rhs.sse) }
299      } else if #[cfg(target_feature="simd128")] {
300        Self { simd: f32x4_ne(self.simd, rhs.simd) }
301      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
302        unsafe {Self { neon: vreinterpretq_f32_u32(vmvnq_u32(vceqq_f32(self.neon, rhs.neon))) }}
303      } else {
304        Self { arr: [
305          if self.arr[0] != rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
306          if self.arr[1] != rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
307          if self.arr[2] != rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
308          if self.arr[3] != rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
309        ]}
310      }
311    }
312  }
313
314  #[inline]
315  fn simd_lt(self, rhs: Self) -> Self::Output {
316    pick! {
317      if #[cfg(target_feature="sse")] {
318        Self { sse: cmp_lt_mask_m128(self.sse, rhs.sse) }
319      } else if #[cfg(target_feature="simd128")] {
320        Self { simd: f32x4_lt(self.simd, rhs.simd) }
321      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
322        unsafe {Self { neon: vreinterpretq_f32_u32(vcltq_f32(self.neon, rhs.neon)) }}
323      } else {
324        Self { arr: [
325          if self.arr[0] < rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
326          if self.arr[1] < rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
327          if self.arr[2] < rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
328          if self.arr[3] < rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
329        ]}
330      }
331    }
332  }
333
334  #[inline]
335  fn simd_gt(self, rhs: Self) -> Self::Output {
336    pick! {
337      if #[cfg(target_feature="sse")] {
338        Self { sse: cmp_gt_mask_m128(self.sse, rhs.sse) }
339      } else if #[cfg(target_feature="simd128")] {
340        Self { simd: f32x4_gt(self.simd, rhs.simd) }
341      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
342        unsafe {Self { neon: vreinterpretq_f32_u32(vcgtq_f32(self.neon, rhs.neon)) }}
343      } else {
344        Self { arr: [
345          if self.arr[0] > rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
346          if self.arr[1] > rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
347          if self.arr[2] > rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
348          if self.arr[3] > rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
349        ]}
350      }
351    }
352  }
353
354  #[inline]
355  fn simd_le(self, rhs: Self) -> Self::Output {
356    pick! {
357      if #[cfg(target_feature="sse")] {
358        Self { sse: cmp_le_mask_m128(self.sse, rhs.sse) }
359      } else if #[cfg(target_feature="simd128")] {
360        Self { simd: f32x4_le(self.simd, rhs.simd) }
361      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
362        unsafe {Self { neon: vreinterpretq_f32_u32(vcleq_f32(self.neon, rhs.neon)) }}
363      } else {
364        Self { arr: [
365          if self.arr[0] <= rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
366          if self.arr[1] <= rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
367          if self.arr[2] <= rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
368          if self.arr[3] <= rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
369        ]}
370      }
371    }
372  }
373
374  #[inline]
375  fn simd_ge(self, rhs: Self) -> Self::Output {
376    pick! {
377      if #[cfg(target_feature="sse")] {
378        Self { sse: cmp_ge_mask_m128(self.sse, rhs.sse) }
379      } else if #[cfg(target_feature="simd128")] {
380        Self { simd: f32x4_ge(self.simd, rhs.simd) }
381      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
382        unsafe {Self { neon: vreinterpretq_f32_u32(vcgeq_f32(self.neon, rhs.neon)) }}
383      } else {
384        Self { arr: [
385          if self.arr[0] >= rhs.arr[0] { f32::from_bits(u32::MAX) } else { 0.0 },
386          if self.arr[1] >= rhs.arr[1] { f32::from_bits(u32::MAX) } else { 0.0 },
387          if self.arr[2] >= rhs.arr[2] { f32::from_bits(u32::MAX) } else { 0.0 },
388          if self.arr[3] >= rhs.arr[3] { f32::from_bits(u32::MAX) } else { 0.0 },
389        ]}
390      }
391    }
392  }
393
394  #[inline]
395  pub fn reduce_add(self) -> f32 {
396    let arr: [f32; 4] = cast(self);
397    arr.iter().sum()
398  }
399
400  #[inline]
401  pub fn reduce_mul(self) -> f32 {
402    let arr: [f32; 4] = cast(self);
403    arr.iter().product()
404  }
405
406  #[inline]
407  pub fn bitselect(self, if_one: Self, if_zero: Self) -> Self {
408    pick! {
409      if #[cfg(target_feature="sse2")] {
410        Self {
411          sse: bitor_m128(
412            bitand_m128(if_one.sse, self.sse),
413            bitandnot_m128(self.sse, if_zero.sse),
414          ),
415        }
416      } else if #[cfg(target_feature="simd128")] {
417        Self { simd: v128_bitselect(if_one.simd, if_zero.simd, self.simd) }
418      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
419        unsafe {Self { neon: vbslq_f32(vreinterpretq_u32_f32(self.neon), if_one.neon, if_zero.neon) }}
420      } else {
421        generic_bit_blend(self, if_one, if_zero)
422      }
423    }
424  }
425
426  #[inline]
427  fn select(self, if_true: Self, if_false: Self) -> Self {
428    pick! {
429      if #[cfg(target_feature="sse4.1")] {
430        Self { sse: blend_varying_m128(if_false.sse, if_true.sse, self.sse) }
431      } else if #[cfg(target_feature="simd128")] {
432        Self { simd: v128_bitselect(if_true.simd, if_false.simd, self.simd) }
433      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
434        unsafe {Self { neon: vbslq_f32(vreinterpretq_u32_f32(self.neon), if_true.neon, if_false.neon) }}
435      } else {
436        generic_bit_blend(self, if_true, if_false)
437      }
438    }
439  }
440
441  #[inline]
442  pub fn to_bitmask(self) -> u32 {
443    pick! {
444      if #[cfg(target_feature="sse")] {
445        move_mask_m128(self.sse) as u32
446      } else if #[cfg(target_feature="simd128")] {
447        u32x4_bitmask(self.simd) as u32
448      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
449        unsafe
450        {
451          // set all to 1 if top bit is set, else 0
452          let masked = vcltq_s32( vreinterpretq_s32_f32(self.neon), vdupq_n_s32(0));
453
454          // select the right bit out of each lane
455          let selectbit : uint32x4_t = core::mem::transmute([1u32, 2, 4, 8]);
456          let r = vandq_u32(masked, selectbit);
457
458          // horizontally add the 16-bit lanes
459          vaddvq_u32(r) as u32
460        }
461      } else {
462        (((self.arr[0].to_bits() as i32) < 0) as u32) |
463        (((self.arr[1].to_bits() as i32) < 0) as u32) << 1 |
464        (((self.arr[2].to_bits() as i32) < 0) as u32) << 2 |
465        (((self.arr[3].to_bits() as i32) < 0) as u32) << 3
466      }
467    }
468  }
469
470  #[inline]
471  pub fn any(self) -> bool {
472    pick! {
473      if #[cfg(target_feature="simd128")] {
474        v128_any_true(self.simd)
475      } else {
476        self.to_bitmask() != 0
477      }
478    }
479  }
480
481  #[inline]
482  pub fn all(self) -> bool {
483    pick! {
484      if #[cfg(target_feature="simd128")] {
485        u32x4_all_true(self.simd)
486      } else {
487        // four lanes
488        self.to_bitmask() == 0b1111
489      }
490    }
491  }
492
493  ///
494  /// Currently this function is only accelerated on `sse`.
495  #[inline]
496  pub fn transpose(data: [f32x4; 4]) -> [f32x4; 4] {
497    pick! {
498      if #[cfg(target_feature="sse")] {
499        let mut e0 = data[0];
500        let mut e1 = data[1];
501        let mut e2 = data[2];
502        let mut e3 = data[3];
503
504        transpose_four_m128(&mut e0.sse, &mut e1.sse, &mut e2.sse, &mut e3.sse);
505
506        [e0, e1, e2, e3]
507      } else if #[cfg(any(all(target_feature="neon",target_arch="aarch64"), target_feature="simd128"))] {
508        let a = data[0].unpack_lo(data[2]);
509        let b = data[1].unpack_lo(data[3]);
510        let c = data[0].unpack_hi(data[2]);
511        let d = data[1].unpack_hi(data[3]);
512
513        [
514          a.unpack_lo(b),
515          a.unpack_hi(b),
516          c.unpack_lo(d),
517          c.unpack_hi(d),
518        ]
519      } else {
520        #[inline(always)]
521        fn transpose_column(data: &[f32x4; 4], index: usize) -> f32x4 {
522          f32x4::new([
523            data[0].as_array()[index],
524            data[1].as_array()[index],
525            data[2].as_array()[index],
526            data[3].as_array()[index],
527          ])
528        }
529
530        [
531          transpose_column(&data, 0),
532          transpose_column(&data, 1),
533          transpose_column(&data, 2),
534          transpose_column(&data, 3),
535        ]
536      }
537    }
538  }
539
540  #[inline]
541  pub fn is_nan(self) -> Self {
542    pick! {
543      if #[cfg(target_feature="sse")] {
544        Self { sse: cmp_unord_mask_m128(self.sse, self.sse) }
545      } else if #[cfg(target_feature="simd128")] {
546        Self { simd: f32x4_ne(self.simd, self.simd) }
547      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
548        unsafe {Self { neon: vreinterpretq_f32_u32(vmvnq_u32(vceqq_f32(self.neon, self.neon))) }}
549      } else {
550        Self { arr: [
551          if self.arr[0].is_nan() { f32::from_bits(u32::MAX) } else { 0.0 },
552          if self.arr[1].is_nan() { f32::from_bits(u32::MAX) } else { 0.0 },
553          if self.arr[2].is_nan() { f32::from_bits(u32::MAX) } else { 0.0 },
554          if self.arr[3].is_nan() { f32::from_bits(u32::MAX) } else { 0.0 },
555        ]}
556      }
557    }
558  }
559
560  #[inline]
561  pub fn is_inf(self) -> Self {
562    let shifted_inf = u32x4::from(0xFF000000);
563    let u: u32x4 = cast(self);
564    let shift_u = u << 1_u64;
565    let out = (shift_u).simd_eq(shifted_inf);
566    cast(out)
567  }
568
569  #[inline]
570  pub fn is_finite(self) -> Self {
571    let shifted_exp_mask = u32x4::from(0xFF000000);
572    let u: u32x4 = cast(self);
573    let shift_u = u << 1_u64;
574    let out = !(shift_u & shifted_exp_mask).simd_eq(shifted_exp_mask);
575    cast(out)
576  }
577
578  #[inline]
579  pub fn is_sign_positive(self) -> Self {
580    const SIGN_MASK: u32x4 = u32x4::splat((-0.0_f32).to_bits());
581
582    let bits = cast::<f32x4, u32x4>(self);
583    let sign = bits & SIGN_MASK;
584    let result = sign.simd_eq(u32x4::ZERO);
585    cast::<u32x4, f32x4>(result)
586  }
587
588  #[inline]
589  pub fn is_sign_negative(self) -> Self {
590    const SIGN_MASK: u32x4 = u32x4::splat((-0.0_f32).to_bits());
591
592    let bits = cast::<f32x4, u32x4>(self);
593    let sign = bits & SIGN_MASK;
594    let result = sign.simd_eq(SIGN_MASK);
595    cast::<u32x4, f32x4>(result)
596  }
597
598  #[inline]
599  pub fn recip(self) -> Self {
600    pick! {
601      if #[cfg(target_feature="sse")] {
602        Self { sse: reciprocal_m128(self.sse) }
603      } else if #[cfg(target_feature="simd128")] {
604        Self { simd: f32x4_div(f32x4_splat(1.0), self.simd) }
605      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
606        unsafe {Self { neon: vdivq_f32(vdupq_n_f32(1.0), self.neon) }}
607      } else {
608        Self { arr: [
609          1.0 / self.arr[0],
610          1.0 / self.arr[1],
611          1.0 / self.arr[2],
612          1.0 / self.arr[3],
613        ]}
614      }
615    }
616  }
617
618  #[inline]
619  pub fn recip_sqrt(self) -> Self {
620    pick! {
621      if #[cfg(target_feature="sse")] {
622        Self { sse: reciprocal_sqrt_m128(self.sse) }
623      } else if #[cfg(target_feature="simd128")] {
624        Self { simd: f32x4_div(f32x4_splat(1.0), f32x4_sqrt(self.simd)) }
625      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
626        unsafe {Self { neon: vdivq_f32(vdupq_n_f32(1.0), vsqrtq_f32(self.neon)) }}
627      } else if #[cfg(feature="std")] {
628        Self { arr: [
629          1.0 / self.arr[0].sqrt(),
630          1.0 / self.arr[1].sqrt(),
631          1.0 / self.arr[2].sqrt(),
632          1.0 / self.arr[3].sqrt(),
633        ]}
634      } else {
635        Self { arr: [
636          1.0 / software_sqrt(self.arr[0] as f64) as f32,
637          1.0 / software_sqrt(self.arr[1] as f64) as f32,
638          1.0 / software_sqrt(self.arr[2] as f64) as f32,
639          1.0 / software_sqrt(self.arr[3] as f64) as f32,
640        ]}
641      }
642    }
643  }
644
645  #[inline]
646  pub fn max(self, rhs: Self) -> Self {
647    pick! {
648      if #[cfg(target_feature="sse")] {
649        // max_m128 seems to do rhs < self ? self : rhs. So if there's any NaN
650        // involved, it chooses rhs, so we need to specifically check rhs for
651        // NaN.
652        rhs.is_nan().select(self, Self { sse: max_m128(self.sse, rhs.sse) })
653      } else if #[cfg(target_feature="simd128")] {
654        // WASM has two max intrinsics:
655        // - max: This propagates NaN, that's the opposite of what we need.
656        // - pmax: This is defined as self < rhs ? rhs : self, which basically
657        //   chooses self if either is NaN.
658        //
659        // pmax is what we want, but we need to specifically check self for NaN.
660        Self {
661          simd: v128_bitselect(
662            rhs.simd,
663            f32x4_pmax(self.simd, rhs.simd),
664            f32x4_ne(self.simd, self.simd), // NaN check
665          )
666        }
667      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
668        unsafe {Self { neon: vmaxnmq_f32(self.neon, rhs.neon) }}
669      } else {
670        Self { arr: [
671          self.arr[0].max(rhs.arr[0]),
672          self.arr[1].max(rhs.arr[1]),
673          self.arr[2].max(rhs.arr[2]),
674          self.arr[3].max(rhs.arr[3]),
675        ]}
676      }
677    }
678  }
679
680  #[inline]
681  pub fn fast_max(self, rhs: Self) -> Self {
682    pick! {
683      if #[cfg(target_feature="sse")] {
684        Self { sse: max_m128(self.sse, rhs.sse) }
685      } else if #[cfg(target_feature="simd128")] {
686        Self {
687          simd: f32x4_pmax(self.simd, rhs.simd),
688        }
689      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
690        unsafe {Self { neon: vmaxq_f32(self.neon, rhs.neon) }}
691      } else {
692        Self { arr: [
693          if self.arr[0] < rhs.arr[0] { rhs.arr[0] } else { self.arr[0] },
694          if self.arr[1] < rhs.arr[1] { rhs.arr[1] } else { self.arr[1] },
695          if self.arr[2] < rhs.arr[2] { rhs.arr[2] } else { self.arr[2] },
696          if self.arr[3] < rhs.arr[3] { rhs.arr[3] } else { self.arr[3] },
697        ]}
698      }
699    }
700  }
701
702  #[inline]
703  pub fn min(self, rhs: Self) -> Self {
704    pick! {
705      if #[cfg(target_feature="sse")] {
706        // min_m128 seems to do self < rhs ? self : rhs. So if there's any NaN
707        // involved, it chooses rhs, so we need to specifically check rhs for
708        // NaN.
709        rhs.is_nan().select(self, Self { sse: min_m128(self.sse, rhs.sse) })
710      } else if #[cfg(target_feature="simd128")] {
711        // WASM has two min intrinsics:
712        // - min: This propagates NaN, that's the opposite of what we need.
713        // - pmin: This is defined as rhs < self ? rhs : self, which basically
714        //   chooses self if either is NaN.
715        //
716        // pmin is what we want, but we need to specifically check self for NaN.
717        Self {
718          simd: v128_bitselect(
719            rhs.simd,
720            f32x4_pmin(self.simd, rhs.simd),
721            f32x4_ne(self.simd, self.simd), // NaN check
722          )
723        }
724      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
725        unsafe {Self { neon: vminnmq_f32(self.neon, rhs.neon) }}
726      } else {
727        Self { arr: [
728          self.arr[0].min(rhs.arr[0]),
729          self.arr[1].min(rhs.arr[1]),
730          self.arr[2].min(rhs.arr[2]),
731          self.arr[3].min(rhs.arr[3]),
732        ]}
733      }
734    }
735  }
736
737  #[inline]
738  pub fn fast_min(self, rhs: Self) -> Self {
739    pick! {
740      if #[cfg(target_feature="sse")] {
741        Self { sse: min_m128(self.sse, rhs.sse) }
742      } else if #[cfg(target_feature="simd128")] {
743        Self {
744          simd: f32x4_pmin(self.simd, rhs.simd),
745        }
746      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
747        unsafe {Self { neon: vminq_f32(self.neon, rhs.neon) }}
748      } else {
749        Self { arr: [
750          if self.arr[0] < rhs.arr[0] { self.arr[0] } else { rhs.arr[0] },
751          if self.arr[1] < rhs.arr[1] { self.arr[1] } else { rhs.arr[1] },
752          if self.arr[2] < rhs.arr[2] { self.arr[2] } else { rhs.arr[2] },
753          if self.arr[3] < rhs.arr[3] { self.arr[3] } else { rhs.arr[3] },
754        ]}
755      }
756    }
757  }
758
759  #[inline]
760  pub fn clamp(self, min: Self, max: Self) -> Self {
761    pick! {
762      if #[cfg(any(
763        target_feature="simd128",
764        all(target_feature="neon",target_arch="aarch64"),
765      ))] {
766        // `fast_clamp` already works.
767        self.fast_clamp(min, max)
768      } else {
769        // This works since all bits set is NaN.
770        self.fast_clamp(min, max) | min.is_nan() | max.is_nan()
771      }
772    }
773  }
774
775  #[inline]
776  pub fn fast_clamp(self, min: Self, max: Self) -> Self {
777    pick! {
778      if #[cfg(target_feature="sse")] {
779        // For both `min_m128` and `max_m128` if any input is NaN, `rhs` gets
780        // chosen. For `self` to be chosen, `self` must be the second argument.
781        Self { sse: max_m128(min.sse, min_m128(max.sse, self.sse)) }
782      } else if #[cfg(target_feature="simd128")] {
783        Self { simd: f32x4_max(f32x4_min(self.simd, max.simd), min.simd) }
784      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
785        unsafe { Self { neon: vmaxq_f32(vminq_f32(self.neon, max.neon), min.neon) } }
786      } else {
787        // The standard library does not have NaN propagating `min` and `max`
788        // functions.
789        let mut result = self;
790        result = result.simd_gt(max).select(max, result);
791        result = result.simd_lt(min).select(min, result);
792        result
793      }
794    }
795  }
796
797  #[inline]
798  pub fn abs(self) -> Self {
799    pick! {
800      if #[cfg(target_feature="simd128")] {
801        Self { simd: f32x4_abs(self.simd) }
802      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
803        unsafe {Self { neon: vabsq_f32(self.neon) }}
804      } else {
805        let non_sign_bits = f32x4::from(f32::from_bits(i32::MAX as u32));
806        self & non_sign_bits
807      }
808    }
809  }
810
811  #[inline]
812  pub fn floor(self) -> Self {
813    pick! {
814      if #[cfg(target_feature="simd128")] {
815        Self { simd: f32x4_floor(self.simd) }
816      } else if #[cfg(target_feature="sse4.1")] {
817        Self { sse: floor_m128(self.sse) }
818      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
819        unsafe {Self { neon: vrndmq_f32(self.neon) }}
820      } else if #[cfg(feature="std")] {
821        let base: [f32; 4] = cast(self);
822        cast(base.map(|val| val.floor()))
823      } else {
824        let base: [f32; 4] = cast(self);
825        let rounded: [f32; 4] = cast(self.round());
826        cast([
827          if base[0] < rounded[0] { rounded[0] - 1.0 } else { rounded[0] },
828          if base[1] < rounded[1] { rounded[1] - 1.0 } else { rounded[1] },
829          if base[2] < rounded[2] { rounded[2] - 1.0 } else { rounded[2] },
830          if base[3] < rounded[3] { rounded[3] - 1.0 } else { rounded[3] },
831        ])
832      }
833    }
834  }
835
836  #[inline]
837  pub fn ceil(self) -> Self {
838    pick! {
839      if #[cfg(target_feature="simd128")] {
840        Self { simd: f32x4_ceil(self.simd) }
841      } else if #[cfg(target_feature="sse4.1")] {
842        Self { sse: ceil_m128(self.sse) }
843      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
844        unsafe {Self { neon: vrndpq_f32(self.neon) }}
845      } else if #[cfg(feature="std")] {
846        let base: [f32; 4] = cast(self);
847        cast(base.map(|val| val.ceil()))
848      } else {
849        let base: [f32; 4] = cast(self);
850        let rounded: [f32; 4] = cast(self.round());
851        cast([
852          if base[0] > rounded[0] { rounded[0] + 1.0 } else { rounded[0] },
853          if base[1] > rounded[1] { rounded[1] + 1.0 } else { rounded[1] },
854          if base[2] > rounded[2] { rounded[2] + 1.0 } else { rounded[2] },
855          if base[3] > rounded[3] { rounded[3] + 1.0 } else { rounded[3] },
856        ])
857      }
858    }
859  }
860
861  #[inline]
862  pub fn round(self) -> Self {
863    pick! {
864      if #[cfg(target_feature="sse4.1")] {
865        const_f32_as_f32x4!(HALF_NEXT_DOWN, 0.5_f32.next_down());
866        const_f32_as_f32x4!(BOUNDS_LIMIT, 8388608.0);
867
868        let self_abs = self.abs();
869
870        let adjusted_self = self_abs + Self::HALF;
871        let result_abs = Self { sse: round_m128::<{round_op!(Zero)}>(adjusted_self.sse) };
872        // The addition breaks for `0.5.next_down()` which incorrectly rounds to
873        // `1.0`. This resets the result back to `0.0`.
874        let result_abs = result_abs & self_abs.simd_ne(HALF_NEXT_DOWN);
875
876        // Large value, infinity and NaN need special handling.
877        let bounds_mask: Self = cast(cmp_lt_mask_i32_m128i(cast(self_abs), cast(BOUNDS_LIMIT)));
878
879        // `abs` keeps the original sign.
880        bounds_mask.abs().bitselect(result_abs, self)
881      } else if #[cfg(target_feature="sse2")] {
882        const_f32_as_f32x4!(HALF_NEXT_DOWN, 0.5_f32.next_down());
883        const_f32_as_f32x4!(BOUNDS_LIMIT, 8388608.0);
884
885        let self_abs = self.abs();
886
887        let adjusted_self = self_abs + Self::HALF;
888        let result_abs = Self {
889          sse: convert_to_m128_from_i32_m128i(truncate_m128_to_m128i(adjusted_self.sse)),
890        };
891        // The addition breaks for `0.5.next_down()` which incorrectly rounds to
892        // `1.0`. This resets the result back to `0.0`.
893        let result_abs = result_abs & self_abs.simd_ne(HALF_NEXT_DOWN);
894
895        // Large value, infinity and NaN need special handling.
896        let bounds_mask: Self = cast(cmp_lt_mask_i32_m128i(cast(self_abs), cast(BOUNDS_LIMIT)));
897
898        // `abs` keeps the original sign.
899        bounds_mask.abs().bitselect(result_abs, self)
900      } else if #[cfg(target_feature="simd128")] {
901        const_f32_as_f32x4!(HALF_NEXT_DOWN, 0.5_f32.next_down());
902        const_f32_as_f32x4!(BOUNDS_LIMIT, 8388608.0);
903
904        let self_abs = self.abs();
905
906        let adjusted_self = self_abs + Self::HALF;
907        let result_abs = Self { simd: f32x4_trunc(adjusted_self.simd) };
908        // The addition breaks for `0.5.next_down()` which incorrectly rounds to
909        // `1.0`. This resets the result back to `0.0`.
910        let result_abs = result_abs & self_abs.simd_ne(HALF_NEXT_DOWN);
911
912        // Large value, infinity and NaN need special handling.
913        let bounds_mask = Self { simd: i32x4_lt(self_abs.simd, BOUNDS_LIMIT.simd) };
914
915        // `abs` keeps the original sign.
916        bounds_mask.abs().bitselect(result_abs, self)
917      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
918        unsafe {Self { neon: vrndaq_f32(self.neon) }}
919      } else {
920        const_f32_as_f32x4!(HALF_NEXT_DOWN, 0.5_f32.next_down());
921        const_f32_as_f32x4!(BOUNDS_LIMIT, 8388608.0);
922
923        let self_abs = self.abs();
924
925        let adjusted_self = (self_abs + Self::HALF).to_array();
926        let result_abs = Self::new([
927          adjusted_self[0] as u32 as f32,
928          adjusted_self[1] as u32 as f32,
929          adjusted_self[2] as u32 as f32,
930          adjusted_self[3] as u32 as f32,
931        ]);
932        // The addition breaks for `0.5.next_down()` which incorrectly rounds to
933        // `1.0`. This resets the result back to `0.0`.
934        let result_abs = result_abs & self_abs.simd_ne(HALF_NEXT_DOWN);
935
936        // Large value, infinity and NaN need special handling.
937        let bounds_mask: Self = cast(cast::<_, i32x4>(self_abs).simd_lt(cast::<_, i32x4>(BOUNDS_LIMIT)));
938
939        // `abs` keeps the original sign.
940        bounds_mask.abs().bitselect(result_abs, self)
941      }
942    }
943  }
944
945  #[inline]
946  pub fn round_int(self) -> i32x4 {
947    pick! {
948      if #[cfg(target_feature="sse2")] {
949        // Based on: https://github.com/v8/v8/blob/210987a552a2bf2a854b0baa9588a5959ff3979d/src/codegen/shared-ia32-x64/macro-assembler-shared-ia32-x64.h#L489-L504
950        let non_nan_mask = self.simd_eq(self);
951        let non_nan = self & non_nan_mask;
952        let flip_to_max: i32x4 = cast(self.simd_ge(Self::splat(2147483648.0)));
953        let cast: i32x4 = cast(convert_to_i32_m128i_from_m128(non_nan.sse));
954        flip_to_max ^ cast
955      } else if #[cfg(target_feature="simd128")] {
956        cast(Self { simd: i32x4_trunc_sat_f32x4(f32x4_nearest(self.simd)) })
957      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
958        cast(unsafe {Self { neon: vreinterpretq_f32_s32(vcvtnq_s32_f32(self.neon)) }})
959      } else {
960        let rounded: [f32; 4] = cast(self.round());
961        cast([
962          rounded[0] as i32,
963          rounded[1] as i32,
964          rounded[2] as i32,
965          rounded[3] as i32,
966        ])
967      }
968    }
969  }
970
971  #[inline]
972  pub fn fast_round_int(self) -> i32x4 {
973    pick! {
974      if #[cfg(target_feature="sse2")] {
975        cast(convert_to_i32_m128i_from_m128(self.sse))
976      } else {
977        self.round_int()
978      }
979    }
980  }
981
982  #[inline]
983  pub fn round_ties_even(self) -> Self {
984    pick! {
985      if #[cfg(target_feature="sse4.1")] {
986        Self { sse: round_m128::<{round_op!(Nearest)}>(self.sse) }
987      } else if #[cfg(target_feature="sse2")] {
988        let mi: m128i = convert_to_i32_m128i_from_m128(self.sse);
989        let f: f32x4 = f32x4 { sse: convert_to_m128_from_i32_m128i(mi) };
990        let i: i32x4 = cast(mi);
991        let mask: f32x4 = cast(i.simd_eq(i32x4::from(0x80000000_u32 as i32)));
992        mask.select(self, f)
993      } else if #[cfg(target_feature="simd128")] {
994        Self { simd: f32x4_nearest(self.simd) }
995      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
996        unsafe {Self { neon: vrndnq_f32(self.neon) }}
997      } else {
998        // Note(Lokathor): This software fallback is probably very slow compared
999        // to having a hardware option available, even just the sse2 version is
1000        // better than this. Oh well.
1001        let to_int = f32x4::from(1.0 / f32::EPSILON);
1002        let u: u32x4 = cast(self);
1003        let e: i32x4 = cast((u >> 23) & u32x4::from(0xff));
1004        let mut y: f32x4;
1005
1006        let no_op_magic = i32x4::from(0x7f + 23);
1007        let no_op_mask: f32x4 = cast(e.simd_gt(no_op_magic) | e.simd_eq(no_op_magic));
1008        let no_op_val: f32x4 = self;
1009
1010        let zero_magic = i32x4::from(0x7f - 1);
1011        let zero_mask: f32x4 = cast(e.simd_lt(zero_magic));
1012        let zero_val: f32x4 = self * f32x4::from(0.0);
1013
1014        let neg_bit: f32x4 = cast(cast::<u32x4, i32x4>(u).simd_lt(i32x4::default()));
1015        let x: f32x4 = neg_bit.select(-self, self);
1016        y = x + to_int - to_int - x;
1017        y = y.simd_gt(f32x4::from(0.5)).select(
1018          y + x - f32x4::from(-1.0),
1019          y.simd_lt(f32x4::from(-0.5)).select(y + x + f32x4::from(1.0), y + x),
1020        );
1021        y = neg_bit.select(-y, y);
1022
1023        no_op_mask.select(no_op_val, zero_mask.select(zero_val, y))
1024      }
1025    }
1026  }
1027
1028  #[inline]
1029  pub fn trunc(self) -> Self {
1030    pick! {
1031      if #[cfg(target_feature="sse4.1")] {
1032        Self { sse: round_m128::<{round_op!(Zero)}>(self.sse) }
1033      } else if #[cfg(target_feature="sse2")] {
1034        // Ported from https://docs.rs/glam/latest/glam/f32/struct.Vec4.html#method.trunc
1035        // Based on https://github.com/microsoft/DirectXMath `XMVectorTruncate`
1036        let result: Self = cast(convert_to_m128_from_i32_m128i(truncate_m128_to_m128i(self.sse)));
1037
1038        // Out of range values are either already round, infinite or NaN.
1039        let bounds_mask: Self = cast(cmp_lt_mask_i32_m128i(
1040            cast(self.abs()),
1041            set_splat_i32_m128i(8388608_f32.to_bits() as i32),
1042        ));
1043
1044        // Reset the sign bit of the mask to preverse the sign of `self`.
1045        bounds_mask.abs().select(result, self)
1046      } else if #[cfg(target_feature="simd128")] {
1047        Self { simd: f32x4_trunc(self.simd) }
1048      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
1049        unsafe { Self { neon: vrndq_f32(self.neon) } }
1050      } else {
1051        let array: [f32; 4] = cast(self);
1052        let result: Self = cast([
1053          array[0] as i32 as f32,
1054          array[1] as i32 as f32,
1055          array[2] as i32 as f32,
1056          array[3] as i32 as f32,
1057        ]);
1058
1059        // Out of range values are either already round, infinite or NaN.
1060        const BOUNDS_LIMIT: i32 = 8388608_f32.to_bits() as i32;
1061        let bounds_mask: Self = cast(cast::<f32x4, i32x4>(self.abs()).simd_lt(i32x4::splat(BOUNDS_LIMIT)));
1062
1063        // Reset the sign bit of the mask to preverse the sign of `self`.
1064        bounds_mask.abs().select(result, self)
1065      }
1066    }
1067  }
1068
1069  #[inline]
1070  pub fn trunc_int(self) -> i32x4 {
1071    pick! {
1072      if #[cfg(target_feature="sse2")] {
1073        // Based on: https://github.com/v8/v8/blob/210987a552a2bf2a854b0baa9588a5959ff3979d/src/codegen/shared-ia32-x64/macro-assembler-shared-ia32-x64.h#L489-L504
1074        let non_nan_mask = self.simd_eq(self);
1075        let non_nan = self & non_nan_mask;
1076        let flip_to_max: i32x4 = cast(self.simd_ge(Self::splat(2147483648.0)));
1077        let cast: i32x4 = cast(truncate_m128_to_m128i(non_nan.sse));
1078        flip_to_max ^ cast
1079      } else if #[cfg(target_feature="simd128")] {
1080        cast(Self { simd: i32x4_trunc_sat_f32x4(self.simd) })
1081      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
1082        cast(unsafe {Self { neon: vreinterpretq_f32_s32(vcvtq_s32_f32(self.neon)) }})
1083      } else {
1084        let n: [f32;4] = cast(self);
1085        cast([
1086          n[0] as i32,
1087          n[1] as i32,
1088          n[2] as i32,
1089          n[3] as i32,
1090        ])
1091      }
1092    }
1093  }
1094
1095  #[inline]
1096  pub fn fast_trunc_int(self) -> i32x4 {
1097    pick! {
1098      if #[cfg(target_feature="sse2")] {
1099        cast(truncate_m128_to_m128i(self.sse))
1100      } else {
1101        self.trunc_int()
1102      }
1103    }
1104  }
1105
1106  ///
1107  /// # Platform-specific behavior (may change in the future)
1108  ///
1109  /// - On `x86`/`x86_64` with FMA: Uses `vfmadd` (single rounding, best
1110  ///   accuracy)
1111  /// - On ARM64 with NEON: Uses `vfmaq_f32` (single rounding, best accuracy)
1112  /// - Without FMA support: Uses `(self * m) + a` (two roundings)
1113  #[inline]
1114  pub fn mul_add(self, a: Self, b: Self) -> Self {
1115    pick! {
1116      if #[cfg(all(target_feature="sse2",target_feature="fma"))] {
1117        Self { sse: fused_mul_add_m128(self.sse, a.sse, b.sse) }
1118      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
1119        unsafe { Self { neon: vfmaq_f32(b.neon, self.neon, a.neon) } }
1120      } else {
1121        (self * a) + b
1122      }
1123    }
1124  }
1125
1126  ///
1127  /// # Platform-specific behavior (may change in the future)
1128  ///
1129  /// - On `x86`/`x86_64` with FMA: Uses `vfmsub` (single rounding, best
1130  ///   accuracy)
1131  /// - On ARM64 with NEON: Uses `vfmaq_f32(-s, self, m)` (single rounding, best
1132  ///   accuracy)
1133  /// - Without FMA support: Uses `(self * m) - s` (two roundings)
1134  #[inline]
1135  pub fn mul_sub(self, a: Self, b: Self) -> Self {
1136    pick! {
1137      if #[cfg(all(target_feature="sse2",target_feature="fma"))] {
1138        Self { sse: fused_mul_sub_m128(self.sse, a.sse, b.sse) }
1139      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
1140        unsafe { Self { neon: vfmaq_f32(vnegq_f32(b.neon), self.neon, a.neon) } }
1141      } else {
1142        (self * a) - b
1143      }
1144    }
1145  }
1146
1147  ///
1148  /// # Platform-specific behavior (may change in the future)
1149  ///
1150  /// - On `x86`/`x86_64` with FMA: Uses `vfnmadd` (single rounding, best
1151  ///   accuracy)
1152  /// - On ARM64 with NEON: Uses `vfmsq_f32` (single rounding, best accuracy)
1153  /// - Without FMA support: Uses `a - (self * m)` (two roundings)
1154  #[inline]
1155  pub fn mul_neg_add(self, a: Self, b: Self) -> Self {
1156    pick! {
1157      if #[cfg(all(target_feature="sse2",target_feature="fma"))] {
1158        Self { sse: fused_mul_neg_add_m128(self.sse, a.sse, b.sse) }
1159      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
1160        unsafe { Self { neon: vfmsq_f32(b.neon, self.neon, a.neon) } }
1161      } else {
1162        b - (self * a)
1163      }
1164    }
1165  }
1166
1167  ///
1168  /// # Platform-specific behavior (may change in the future)
1169  ///
1170  /// - On `x86`/`x86_64` with FMA: Uses `vfnmsub` (single rounding, best
1171  ///   accuracy)
1172  /// - On ARM64 with NEON: Uses `-(vfmaq_f32(s, self, m))` (single rounding,
1173  ///   best accuracy)
1174  /// - Without FMA support: Uses `-(self * m) - s` (two roundings)
1175  #[inline]
1176  pub fn mul_neg_sub(self, a: Self, b: Self) -> Self {
1177    pick! {
1178      if #[cfg(all(target_feature="sse2",target_feature="fma"))] {
1179        Self { sse: fused_mul_neg_sub_m128(self.sse, a.sse, b.sse) }
1180      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
1181        unsafe { Self { neon: vnegq_f32(vfmaq_f32(b.neon, self.neon, a.neon)) } }
1182      } else {
1183        -(self * a) - b
1184      }
1185    }
1186  }
1187
1188  #[inline]
1189  pub fn powf_simd(self, n: Self) -> Self {
1190    const_f32_as_f32x4!(ln2f_hi, 0.693359375);
1191    const_f32_as_f32x4!(ln2f_lo, -2.12194440e-4);
1192    const_f32_as_f32x4!(P0logf, 3.3333331174E-1);
1193    const_f32_as_f32x4!(P1logf, -2.4999993993E-1);
1194    const_f32_as_f32x4!(P2logf, 2.0000714765E-1);
1195    const_f32_as_f32x4!(P3logf, -1.6668057665E-1);
1196    const_f32_as_f32x4!(P4logf, 1.4249322787E-1);
1197    const_f32_as_f32x4!(P5logf, -1.2420140846E-1);
1198    const_f32_as_f32x4!(P6logf, 1.1676998740E-1);
1199    const_f32_as_f32x4!(P7logf, -1.1514610310E-1);
1200    const_f32_as_f32x4!(P8logf, 7.0376836292E-2);
1201
1202    const_f32_as_f32x4!(p2expf, 1.0 / 2.0); // coefficients for Taylor expansion of exp
1203    const_f32_as_f32x4!(p3expf, 1.0 / 6.0);
1204    const_f32_as_f32x4!(p4expf, 1.0 / 24.0);
1205    const_f32_as_f32x4!(p5expf, 1.0 / 120.0);
1206    const_f32_as_f32x4!(p6expf, 1.0 / 720.0);
1207    const_f32_as_f32x4!(p7expf, 1.0 / 5040.0);
1208
1209    let x1 = self.abs();
1210    let x = x1.fraction_2();
1211
1212    let mask = x.simd_gt(f32x4::SQRT_2 * f32x4::HALF);
1213    let x = (!mask).select(x + x, x);
1214
1215    let x = x - f32x4::ONE;
1216    let x2 = x * x;
1217    let lg1 = polynomial_8!(
1218      x, P0logf, P1logf, P2logf, P3logf, P4logf, P5logf, P6logf, P7logf, P8logf
1219    );
1220    let lg1 = lg1 * x2 * x;
1221
1222    let ef = x1.exponent();
1223    let ef = mask.select(ef + f32x4::ONE, ef);
1224
1225    let e1 = (ef * n).round_ties_even();
1226    let yr = ef.mul_sub(n, e1);
1227
1228    let lg = f32x4::HALF.mul_neg_add(x2, x) + lg1;
1229    let x2_err = (f32x4::HALF * x).mul_sub(x, f32x4::HALF * x2);
1230    let lg_err = f32x4::HALF.mul_add(x2, lg - x) - lg1;
1231
1232    let e2 = (lg * n * f32x4::LOG2_E).round_ties_even();
1233    let v = lg.mul_sub(n, e2 * ln2f_hi);
1234    let v = e2.mul_neg_add(ln2f_lo, v);
1235    let v = v - (lg_err + x2_err).mul_sub(n, yr * f32x4::LN_2);
1236
1237    let x = v;
1238    let e3 = (x * f32x4::LOG2_E).round_ties_even();
1239    let x = e3.mul_neg_add(f32x4::LN_2, x);
1240    let x2 = x * x;
1241    let z = x2.mul_add(
1242      polynomial_5!(x, p2expf, p3expf, p4expf, p5expf, p6expf, p7expf),
1243      x + f32x4::ONE,
1244    );
1245
1246    let ee = e1 + e2 + e3;
1247    let ei = cast::<_, i32x4>(ee.round_int());
1248    let ej = cast::<_, i32x4>(ei + (cast::<_, i32x4>(z) >> 23));
1249
1250    let overflow = cast::<_, f32x4>(ej.simd_gt(i32x4::splat(0x0FF)))
1251      | (ee.simd_gt(f32x4::splat(300.0)));
1252    let underflow = cast::<_, f32x4>(ej.simd_lt(i32x4::splat(0x000)))
1253      | (ee.simd_lt(f32x4::splat(-300.0)));
1254
1255    // Add exponent by integer addition
1256    let z = cast::<_, f32x4>(cast::<_, i32x4>(z) + (ei << 23));
1257
1258    // Check for overflow/underflow
1259    let z = if (overflow | underflow).any() {
1260      let z = underflow.select(f32x4::ZERO, z);
1261      overflow.select(Self::infinity(), z)
1262    } else {
1263      z
1264    };
1265
1266    // Check for self == 0
1267    let x_zero = self.is_zero_or_subnormal();
1268    let z = x_zero.select(
1269      n.simd_lt(f32x4::ZERO).select(
1270        Self::infinity(),
1271        n.simd_eq(f32x4::ZERO).select(f32x4::ONE, f32x4::ZERO),
1272      ),
1273      z,
1274    );
1275
1276    let x_sign = self.is_sign_negative();
1277    let z = if x_sign.any() {
1278      // Y into an integer
1279      let yi = n.simd_eq(n.round_ties_even());
1280      // Is y odd? If yes flip the sign of the result.
1281      let y_odd = cast::<i32x4, f32x4>(n.round_int() << 31);
1282
1283      let z1 = yi
1284        .select(z | y_odd, self.simd_eq(Self::ZERO).select(z, Self::nan_pow()));
1285      x_sign.select(z1, z)
1286    } else {
1287      z
1288    };
1289
1290    let x_finite = self.is_finite();
1291    let y_finite = n.is_finite();
1292    let e_finite = ee.is_finite();
1293    if (x_finite & y_finite & (e_finite | x_zero)).all() {
1294      return z;
1295    }
1296
1297    (self.is_nan() | n.is_nan()).select(self + n, z)
1298  }
1299
1300  #[inline]
1301  pub fn sqrt(self) -> Self {
1302    pick! {
1303      if #[cfg(target_feature="sse")] {
1304        Self { sse: sqrt_m128(self.sse) }
1305      } else if #[cfg(target_feature="simd128")] {
1306        Self { simd: f32x4_sqrt(self.simd) }
1307      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
1308        unsafe {Self { neon: vsqrtq_f32(self.neon) }}
1309      } else if #[cfg(feature="std")] {
1310        Self { arr: [
1311          self.arr[0].sqrt(),
1312          self.arr[1].sqrt(),
1313          self.arr[2].sqrt(),
1314          self.arr[3].sqrt(),
1315        ]}
1316      } else {
1317        Self { arr: [
1318          software_sqrt(self.arr[0] as f64) as f32,
1319          software_sqrt(self.arr[1] as f64) as f32,
1320          software_sqrt(self.arr[2] as f64) as f32,
1321          software_sqrt(self.arr[3] as f64) as f32,
1322        ]}
1323      }
1324    }
1325  }
1326
1327  /// Calculate the exponent of a packed `f32x4`
1328  #[inline]
1329  pub fn exp(self) -> Self {
1330    const_f32_as_f32x4!(P0, 1.0 / 2.0);
1331    const_f32_as_f32x4!(P1, 1.0 / 6.0);
1332    const_f32_as_f32x4!(P2, 1.0 / 24.0);
1333    const_f32_as_f32x4!(P3, 1.0 / 120.0);
1334    const_f32_as_f32x4!(P4, 1.0 / 720.0);
1335    const_f32_as_f32x4!(P5, 1.0 / 5040.0);
1336    // LN2D_HI/LO: double-double decomposition of ln(2) for exp range reduction,
1337    // following the approach from fdlibm's e_exp.c (Sun Microsystems,
1338    // https://www.netlib.org/fdlibm/). The f32 split uses f32-precision constants
1339    // (0.693359375, -2.12194440e-4) summing to ln(2) with single-precision
1340    // accuracy; the f64 variants use a full f64 double-double
1341    // decomposition.
1342    const_f32_as_f32x4!(LN2D_HI, 0.693359375);
1343    const_f32_as_f32x4!(LN2D_LO, -2.12194440e-4);
1344    // max_x = ln(f32::MAX) ≈ 88.7229, max_r = 127 (IEEE max normal exponent)
1345    // min_x = -149.5 ln(2) ≈ -103.63: min r for vm_pow2n subnormal
1346    let max_x = f32x4::from(88.723);
1347    let min_x = f32x4::from(-103.63);
1348    // x < min_x: e^x underflows to 0 -- skip the entire pipeline
1349    let finite = self.is_finite();
1350    let neg_underflow = self.simd_lt(min_x) & finite;
1351    if neg_underflow.all() {
1352      return Self::ZERO;
1353    }
1354    let max_r = f32x4::from(127.0);
1355    let r = (self * Self::LOG2_E).round_ties_even();
1356    let big = r.simd_gt(max_r);
1357    let r_safe = big.select(max_r, r);
1358    let excess = r - max_r;
1359    let excess = big.select(excess, Self::ZERO);
1360    let scale = Self::vm_pow2n(excess);
1361    let x = r.mul_neg_add(LN2D_HI, self);
1362    let x = r.mul_neg_add(LN2D_LO, x);
1363    let z = polynomial_5!(x, P0, P1, P2, P3, P4, P5);
1364    let x2 = x * x;
1365    let z = z.mul_add(x2, x);
1366    let n2 = Self::vm_pow2n(r_safe);
1367    let z = (z + Self::ONE) * scale * n2;
1368    let nan_mask = self.is_nan();
1369    let mut result = nan_mask.select(Self::nan_pow(), z);
1370    let pos_overflow = self.simd_gt(max_x) & finite;
1371    result = pos_overflow.select(Self::infinity(), result);
1372    result = neg_underflow.select(Self::ZERO, result);
1373    let pos_inf = !finite & !self.is_sign_negative() & !nan_mask;
1374    result = pos_inf.select(Self::infinity(), result);
1375    let neg_inf = !finite & self.is_sign_negative() & !nan_mask;
1376    result = neg_inf.select(Self::ZERO, result);
1377    result
1378  }
1379
1380  #[inline]
1381  pub fn exp2(self) -> Self {
1382    const_f32_as_f32x4!(P2, 1.0 / 2.0);
1383    const_f32_as_f32x4!(P3, 1.0 / 6.0);
1384    const_f32_as_f32x4!(P4, 1.0 / 24.0);
1385    const_f32_as_f32x4!(P5, 1.0 / 120.0);
1386    const_f32_as_f32x4!(P6, 1.0 / 720.0);
1387    const_f32_as_f32x4!(P7, 1.0 / 5040.0);
1388
1389    // max_x = log2(f32::MAX) ≈ 127.99999
1390    // min_x = log2(f32::MIN_POSITIVE) - 23 ≈ -126 - 23 = -149
1391    let max_x = f32x4::from(127.99999);
1392    let min_x = f32x4::from(-149.5);
1393    let finite = self.is_finite();
1394    let neg_underflow = self.simd_lt(min_x) & finite;
1395    if neg_underflow.all() {
1396      return Self::ZERO;
1397    }
1398
1399    let round = self.round_ties_even();
1400    let max_r = f32x4::from(127.0);
1401    let big = round.simd_gt(max_r);
1402    let r_safe = big.select(max_r, round);
1403    let excess = round - max_r;
1404    let excess = big.select(excess, Self::ZERO);
1405    let scale = Self::vm_pow2n(excess);
1406
1407    let fract = (self - round) * Self::LN_2;
1408    let fract_partial_exp2 = polynomial_5!(fract, P2, P3, P4, P5, P6, P7);
1409    let fract2 = fract * fract;
1410    let fract_exp2 = fract_partial_exp2.mul_add(fract2, fract) + Self::ONE;
1411
1412    let n2 = Self::vm_pow2n(r_safe);
1413    let result = fract_exp2 * scale * n2;
1414
1415    let nan_mask = self.is_nan();
1416    let mut result = nan_mask.select(Self::nan_pow(), result);
1417    let pos_overflow = self.simd_gt(max_x) & finite;
1418    result = pos_overflow.select(Self::infinity(), result);
1419    result = neg_underflow.select(Self::ZERO, result);
1420    let pos_inf = !finite & !self.is_sign_negative() & !nan_mask;
1421    result = pos_inf.select(Self::infinity(), result);
1422    let neg_inf = !finite & self.is_sign_negative() & !nan_mask;
1423    result = neg_inf.select(Self::ZERO, result);
1424    result
1425  }
1426
1427  #[inline]
1428  pub fn ln(self) -> Self {
1429    const_f32_as_f32x4!(HALF, 0.5);
1430    const_f32_as_f32x4!(P0, 3.3333331174E-1);
1431    const_f32_as_f32x4!(P1, -2.4999993993E-1);
1432    const_f32_as_f32x4!(P2, 2.0000714765E-1);
1433    const_f32_as_f32x4!(P3, -1.6668057665E-1);
1434    const_f32_as_f32x4!(P4, 1.4249322787E-1);
1435    const_f32_as_f32x4!(P5, -1.2420140846E-1);
1436    const_f32_as_f32x4!(P6, 1.1676998740E-1);
1437    const_f32_as_f32x4!(P7, -1.1514610310E-1);
1438    const_f32_as_f32x4!(P8, 7.0376836292E-2);
1439    const_f32_as_f32x4!(LN2F_HI, 0.693359375);
1440    const_f32_as_f32x4!(LN2F_LO, -2.12194440e-4);
1441    const_f32_as_f32x4!(VM_SMALLEST_NORMAL, 1.17549435E-38);
1442
1443    let x1 = self;
1444    let x = Self::fraction_2(x1);
1445    let e = Self::exponent(x1);
1446    let mask = x.simd_gt(Self::SQRT_2 * HALF);
1447    let x = (!mask).select(x + x, x);
1448    let fe = mask.select(e + Self::ONE, e);
1449    let x = x - Self::ONE;
1450    let res = polynomial_8!(x, P0, P1, P2, P3, P4, P5, P6, P7, P8);
1451    let x2 = x * x;
1452    let res = x2 * x * res;
1453    let res = fe.mul_add(LN2F_LO, res);
1454    let res = res + x2.mul_neg_add(HALF, x);
1455    let res = fe.mul_add(LN2F_HI, res);
1456    let overflow = !self.is_finite();
1457    let underflow = x1.simd_lt(VM_SMALLEST_NORMAL);
1458    let mask = overflow | underflow;
1459    if !mask.any() {
1460      res
1461    } else {
1462      let is_zero = self.is_zero_or_subnormal();
1463      let res = underflow.select(Self::nan_log(), res);
1464      // Note: is_zero_or_subnormal() lumps subnormals (exponent==0) with zero.
1465      // Both get -Inf here. True subnormal inputs (~1.4e-45..1.175e-38) should
1466      // produce a finite negative result, but are vanishingly rare in
1467      // practice.
1468      let res = is_zero.select(-Self::infinity(), res);
1469      let res = overflow.select(self, res);
1470      // This must come *after* overflow.blend to overwrite ln(-∞) = -∞ to NaN
1471      let res = (!self.is_finite() & self.is_sign_negative())
1472        .select(Self::nan_log(), res);
1473      res
1474    }
1475  }
1476
1477  #[inline]
1478  pub fn cbrt(self) -> Self {
1479    let a = self.abs();
1480    let zero = a.simd_eq(Self::ZERO);
1481    if zero.all() {
1482      return self; // preserves -0.0
1483    }
1484    let inf = a.is_inf();
1485    let nan = self.is_nan();
1486
1487    let tiny = a.simd_lt(Self::from(f32::MIN_POSITIVE));
1488    let a_work = tiny.select(a * Self::from(16777216.0), a);
1489
1490    let e = Self::exponent(a_work) + Self::ONE;
1491    let d = Self::fraction_2(a_work);
1492
1493    // C0..C5 from SLEEF's minimax polynomial for 1/cbrt(d) on [0.5, 1.0)
1494    // Naoki Shibata et al., "SLEEF: A Portable Vectorized Library of C99
1495    // Mathematical Functions", https://sleef.org / https://github.com/shibatch/sleef
1496    // Licensed under the Boost Software License 1.0.
1497    // These are the f32-precision coefficients; our f64 variants use the f64
1498    // set.
1499    const_f32_as_f32x4!(C0, 2.2241257);
1500    const_f32_as_f32x4!(C1, -3.8095417);
1501    const_f32_as_f32x4!(C2, 5.8982625);
1502    const_f32_as_f32x4!(C3, -5.532182);
1503    const_f32_as_f32x4!(C4, 2.8208892);
1504    const_f32_as_f32x4!(C5, -0.60156447);
1505    let mut x = polynomial_5!(d, C0, C1, C2, C3, C4, C5);
1506
1507    let x2 = x * x;
1508    let x4 = x2 * x2;
1509    x = x - d.mul_add(x4, -x) * Self::from(1.0 / 3.0);
1510    // cbrt(d) = d * x² with refinement
1511    let mut y = (d * x) * x;
1512    let yx = y * x;
1513    let t = Self::from(2.0 / 3.0);
1514    y = y - t * y * (yx - Self::ONE);
1515
1516    // Scale by 2^(e/3)
1517    let three = Self::from(3.0);
1518    let two = Self::from(2.0);
1519    let neg = e.simd_lt(Self::ZERO);
1520    let e_adj = neg.select(e - two, e);
1521    let k = (e_adj / three).trunc();
1522    let r = e - three * k;
1523    const_f32_as_f32x4!(CBRT2, 1.259921);
1524    const_f32_as_f32x4!(CBRT4, 1.587401);
1525    y = r.simd_eq(Self::ONE).select(y * CBRT2, y);
1526    y = r.simd_eq(two).select(y * CBRT4, y);
1527    y *= Self::vm_pow2n(k);
1528    y = tiny.select(y / Self::from(256.0_f32), y);
1529
1530    let result = y.flip_signs(self);
1531    let result = nan.select(self, result);
1532    let result = zero.select(self, result);
1533    let result = inf.select(self, result);
1534    result
1535  }
1536
1537  #[inline]
1538  pub fn asin(self) -> Self {
1539    // Based on the Agner Fog "vector class library":
1540    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
1541    const_f32_as_f32x4!(P4asinf, 4.2163199048E-2);
1542    const_f32_as_f32x4!(P3asinf, 2.4181311049E-2);
1543    const_f32_as_f32x4!(P2asinf, 4.5470025998E-2);
1544    const_f32_as_f32x4!(P1asinf, 7.4953002686E-2);
1545    const_f32_as_f32x4!(P0asinf, 1.6666752422E-1);
1546
1547    let xa = self.abs();
1548    let big = xa.simd_ge(f32x4::splat(0.5));
1549
1550    let x1 = f32x4::splat(0.5) * (f32x4::ONE - xa);
1551    let x2 = xa * xa;
1552    let x3 = big.select(x1, x2);
1553
1554    let xb = x1.sqrt();
1555
1556    let x4 = big.select(xb, xa);
1557
1558    let z = polynomial_4!(x3, P0asinf, P1asinf, P2asinf, P3asinf, P4asinf);
1559    let z = z.mul_add(x3 * x4, x4);
1560
1561    let z1 = z + z;
1562
1563    // asin
1564    let z3 = f32x4::FRAC_PI_2 - z1;
1565    let asin = big.select(z3, z);
1566    let asin = asin.flip_signs(self);
1567
1568    asin
1569  }
1570
1571  #[inline]
1572  pub fn acos(self) -> Self {
1573    // Based on the Agner Fog "vector class library":
1574    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
1575    const_f32_as_f32x4!(P4asinf, 4.2163199048E-2);
1576    const_f32_as_f32x4!(P3asinf, 2.4181311049E-2);
1577    const_f32_as_f32x4!(P2asinf, 4.5470025998E-2);
1578    const_f32_as_f32x4!(P1asinf, 7.4953002686E-2);
1579    const_f32_as_f32x4!(P0asinf, 1.6666752422E-1);
1580
1581    let xa = self.abs();
1582    let big = xa.simd_ge(f32x4::splat(0.5));
1583
1584    let x1 = f32x4::splat(0.5) * (f32x4::ONE - xa);
1585    let x2 = xa * xa;
1586    let x3 = big.select(x1, x2);
1587
1588    let xb = x1.sqrt();
1589
1590    let x4 = big.select(xb, xa);
1591
1592    let z = polynomial_4!(x3, P0asinf, P1asinf, P2asinf, P3asinf, P4asinf);
1593    let z = z.mul_add(x3 * x4, x4);
1594
1595    let z1 = z + z;
1596
1597    // acos
1598    let z3 = self.simd_lt(f32x4::ZERO).select(f32x4::PI - z1, z1);
1599    let z4 = f32x4::FRAC_PI_2 - z.flip_signs(self);
1600    let acos = big.select(z3, z4);
1601
1602    acos
1603  }
1604
1605  #[inline]
1606  pub fn atan(self) -> Self {
1607    // Based on the Agner Fog "vector class library":
1608    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
1609    const_f32_as_f32x4!(P3atanf, 8.05374449538E-2);
1610    const_f32_as_f32x4!(P2atanf, -1.38776856032E-1);
1611    const_f32_as_f32x4!(P1atanf, 1.99777106478E-1);
1612    const_f32_as_f32x4!(P0atanf, -3.33329491539E-1);
1613
1614    let t = self.abs();
1615
1616    // small:  z = t / 1.0;
1617    // medium: z = (t-1.0) / (t+1.0);
1618    // big:    z = -1.0 / t;
1619    let notsmal = t.simd_ge(Self::SQRT_2 - Self::ONE);
1620    let notbig = t.simd_le(Self::SQRT_2 + Self::ONE);
1621
1622    let mut s = notbig.select(Self::FRAC_PI_4, Self::FRAC_PI_2);
1623    s = notsmal & s;
1624
1625    let mut a = notbig & t;
1626    a = notsmal.select(a - Self::ONE, a);
1627    let mut b = notbig & Self::ONE;
1628    b = notsmal.select(b + t, b);
1629    let z = a / b;
1630
1631    let zz = z * z;
1632
1633    // Taylor expansion
1634    let mut re = polynomial_3!(zz, P0atanf, P1atanf, P2atanf, P3atanf);
1635    re = re.mul_add(zz * z, z) + s;
1636
1637    // get sign bit
1638    re = (self.is_sign_negative()).select(-re, re);
1639
1640    re
1641  }
1642
1643  #[inline]
1644  pub fn atan2(self, x: Self) -> Self {
1645    // Based on the Agner Fog "vector class library":
1646    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
1647    const_f32_as_f32x4!(P3atanf, 8.05374449538E-2);
1648    const_f32_as_f32x4!(P2atanf, -1.38776856032E-1);
1649    const_f32_as_f32x4!(P1atanf, 1.99777106478E-1);
1650    const_f32_as_f32x4!(P0atanf, -3.33329491539E-1);
1651
1652    let y = self;
1653
1654    // move in first octant
1655    let x1 = x.abs();
1656    let y1 = y.abs();
1657    let swapxy = y1.simd_gt(x1);
1658    // swap x and y if y1 > x1
1659    let mut x2 = swapxy.select(y1, x1);
1660    let mut y2 = swapxy.select(x1, y1);
1661
1662    // check for special case: x and y are both +/- INF
1663    let both_infinite = x.is_inf() & y.is_inf();
1664    if both_infinite.any() {
1665      let minus_one = -Self::ONE;
1666      x2 = both_infinite.select(x2 & minus_one, x2);
1667      y2 = both_infinite.select(y2 & minus_one, y2);
1668    }
1669
1670    // x = y = 0 will produce NAN. No problem, fixed below
1671    let t = y2 / x2;
1672
1673    // small:  z = t / 1.0;
1674    // medium: z = (t-1.0) / (t+1.0);
1675    let notsmal = t.simd_ge(Self::SQRT_2 - Self::ONE);
1676
1677    let a = notsmal.select(t - Self::ONE, t);
1678    let b = notsmal.select(t + Self::ONE, Self::ONE);
1679    let s = notsmal & Self::FRAC_PI_4;
1680    let z = a / b;
1681
1682    let zz = z * z;
1683
1684    // Taylor expansion
1685    let mut re = polynomial_3!(zz, P0atanf, P1atanf, P2atanf, P3atanf);
1686    re = re.mul_add(zz * z, z) + s;
1687
1688    // move back in place
1689    re = swapxy.select(Self::FRAC_PI_2 - re, re);
1690    re = ((x | y).simd_eq(Self::ZERO)).select(Self::ZERO, re);
1691    re = (x.is_sign_negative()).select(Self::PI - re, re);
1692
1693    // get sign bit
1694    re = (y.is_sign_negative()).select(-re, re);
1695
1696    re
1697  }
1698
1699  #[inline]
1700  pub fn sin_cos(self) -> (Self, Self) {
1701    // Based on the Agner Fog "vector class library":
1702    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
1703
1704    const_f32_as_f32x4!(DP1F, 0.78515625_f32 * 2.0);
1705    const_f32_as_f32x4!(DP2F, 2.4187564849853515625E-4_f32 * 2.0);
1706    const_f32_as_f32x4!(DP3F, 3.77489497744594108E-8_f32 * 2.0);
1707
1708    const_f32_as_f32x4!(P0sinf, -1.6666654611E-1);
1709    const_f32_as_f32x4!(P1sinf, 8.3321608736E-3);
1710    const_f32_as_f32x4!(P2sinf, -1.9515295891E-4);
1711
1712    const_f32_as_f32x4!(P0cosf, 4.166664568298827E-2);
1713    const_f32_as_f32x4!(P1cosf, -1.388731625493765E-3);
1714    const_f32_as_f32x4!(P2cosf, 2.443315711809948E-5);
1715
1716    const_f32_as_f32x4!(TWO_OVER_PI, 2.0 / core::f32::consts::PI);
1717
1718    let xa = self.abs();
1719
1720    // Find quadrant
1721    let y = (xa * TWO_OVER_PI).round_ties_even();
1722    let q: i32x4 = y.round_int();
1723
1724    let x = y.mul_neg_add(DP3F, y.mul_neg_add(DP2F, y.mul_neg_add(DP1F, xa)));
1725
1726    let x2 = x * x;
1727    let mut s = polynomial_2!(x2, P0sinf, P1sinf, P2sinf) * (x * x2) + x;
1728    let mut c = polynomial_2!(x2, P0cosf, P1cosf, P2cosf) * (x2 * x2)
1729      + f32x4::from(0.5).mul_neg_add(x2, f32x4::from(1.0));
1730
1731    let swap = !(q & i32x4::from(1)).simd_eq(i32x4::from(0));
1732
1733    let mut overflow: f32x4 = cast(q.simd_gt(i32x4::from(0x2000000)));
1734    overflow &= xa.is_finite();
1735    s = overflow.select(f32x4::from(0.0), s);
1736    c = overflow.select(f32x4::from(1.0), c);
1737
1738    // calc sin
1739    let mut sin1 = cast::<_, f32x4>(swap).select(c, s);
1740    let sign_sin: i32x4 = (q << 30) ^ cast::<_, i32x4>(self);
1741    sin1 = sin1.flip_signs(cast(sign_sin));
1742
1743    // calc cos
1744    let mut cos1 = cast::<_, f32x4>(swap).select(s, c);
1745    let sign_cos: i32x4 = ((q + i32x4::from(1)) & i32x4::from(2)) << 30;
1746    cos1 ^= cast::<_, f32x4>(sign_cos);
1747
1748    // IEEE 754: sin/cos(±∞) = NaN, sin/cos(NaN) = NaN
1749    let finite = self.is_finite();
1750    let nan = Self::splat(f32::NAN);
1751    let sin_final = finite.select(sin1, nan);
1752    let cos_final = finite.select(cos1, nan);
1753
1754    (sin_final, cos_final)
1755  }
1756
1757  #[inline]
1758  pub fn asin_acos(self) -> (Self, Self) {
1759    // Based on the Agner Fog "vector class library":
1760    // https://github.com/vectorclass/version2/blob/master/vectormath_trig.h
1761    const_f32_as_f32x4!(P4asinf, 4.2163199048E-2);
1762    const_f32_as_f32x4!(P3asinf, 2.4181311049E-2);
1763    const_f32_as_f32x4!(P2asinf, 4.5470025998E-2);
1764    const_f32_as_f32x4!(P1asinf, 7.4953002686E-2);
1765    const_f32_as_f32x4!(P0asinf, 1.6666752422E-1);
1766
1767    let xa = self.abs();
1768    let big = xa.simd_ge(f32x4::splat(0.5));
1769
1770    let x1 = f32x4::splat(0.5) * (f32x4::ONE - xa);
1771    let x2 = xa * xa;
1772    let x3 = big.select(x1, x2);
1773
1774    let xb = x1.sqrt();
1775
1776    let x4 = big.select(xb, xa);
1777
1778    let z = polynomial_4!(x3, P0asinf, P1asinf, P2asinf, P3asinf, P4asinf);
1779    let z = z.mul_add(x3 * x4, x4);
1780
1781    let z1 = z + z;
1782
1783    // acos
1784    let z3 = self.simd_lt(f32x4::ZERO).select(f32x4::PI - z1, z1);
1785    let z4 = f32x4::FRAC_PI_2 - z.flip_signs(self);
1786    let acos = big.select(z3, z4);
1787
1788    // asin
1789    let z3 = f32x4::FRAC_PI_2 - z1;
1790    let asin = big.select(z3, z);
1791    let asin = asin.flip_signs(self);
1792
1793    (asin, acos)
1794  }
1795
1796  #[inline]
1797  pub fn exp_m1(self) -> Self {
1798    // x < -17.329: e^x < 2⁻²⁵, exp_m1(x) = -1.0 exactly (mantissa exhaustion)
1799    // IEEE simd_lt returns false for NaN, so NaN lanes can't reach here.
1800    // -inf is < -17.329, and exp_m1(-inf) = -1.0, also correct.
1801    if self.simd_lt(f32x4::from(-17.329)).all() {
1802      return f32x4::from(-1.0);
1803    }
1804    const_f32_as_f32x4!(P0, 1.0 / 2.0);
1805    const_f32_as_f32x4!(P1, 1.0 / 6.0);
1806    const_f32_as_f32x4!(P2, 1.0 / 24.0);
1807    const_f32_as_f32x4!(P3, 1.0 / 120.0);
1808    const_f32_as_f32x4!(P4, 1.0 / 720.0);
1809    const_f32_as_f32x4!(P5, 1.0 / 5040.0);
1810    // LN2D_HI/LO: double-double decomposition of ln(2) for exp range reduction,
1811    // following the approach from fdlibm's e_exp.c (Sun Microsystems,
1812    // https://www.netlib.org/fdlibm/). The f32 split uses f32-precision constants
1813    // (0.693359375, -2.12194440e-4) summing to ln(2) with single-precision
1814    // accuracy; the f64 variants use a full f64 double-double
1815    // decomposition.
1816    const_f32_as_f32x4!(LN2D_HI, 0.693359375);
1817    const_f32_as_f32x4!(LN2D_LO, -2.12194440e-4);
1818    // max_x = ln(f32::MAX) ≈ 88.7229, max_r = 127 (IEEE max normal exponent)
1819    // min_x = -149.5 ln(2) ≈ -103.63: min r for vm_pow2n subnormal
1820    let max_x = f32x4::from(88.723);
1821    let min_x = f32x4::from(-103.63);
1822    let max_r = f32x4::from(127.0);
1823    let r = (self * Self::LOG2_E).round_ties_even();
1824    let big = r.simd_gt(max_r);
1825    let r_safe = big.select(max_r, r);
1826    let excess = r - max_r;
1827    let excess = big.select(excess, Self::ZERO);
1828    let scale = Self::vm_pow2n(excess);
1829    let x = r.mul_neg_add(LN2D_HI, self);
1830    let x = r.mul_neg_add(LN2D_LO, x);
1831    let z = polynomial_5!(x, P0, P1, P2, P3, P4, P5);
1832    let x2 = x * x;
1833    let z = z.mul_add(x2, x);
1834    let n2 = Self::vm_pow2n(r_safe);
1835    let exp_val = (z + Self::ONE) * scale * n2;
1836    let r_is_zero = r.simd_eq(Self::ZERO);
1837    let z = r_is_zero.select(z, exp_val - Self::ONE);
1838    let nan_mask = self.is_nan();
1839    let finite = self.is_finite();
1840    let mut result = nan_mask.select(Self::nan_pow(), z);
1841    let pos_overflow = self.simd_gt(max_x) & finite;
1842    result = pos_overflow.select(Self::infinity(), result);
1843    let neg_underflow = self.simd_lt(min_x) & finite;
1844    result = neg_underflow.select(-Self::ONE, result);
1845    let pos_inf = !finite & !self.is_sign_negative() & !nan_mask;
1846    result = pos_inf.select(Self::infinity(), result);
1847    let neg_inf = !finite & self.is_sign_negative() & !nan_mask;
1848    result = neg_inf.select(-Self::ONE, result);
1849    let is_zero = self.simd_eq(Self::ZERO);
1850    result = is_zero.select(self, result);
1851    result
1852  }
1853
1854  #[inline]
1855  pub fn ln_1p(self) -> Self {
1856    // Based on the identity ln(1+x) = x·ln(1+x)/((1+x)-1), i.e. x·ln(u)/(u-1)
1857    // where u = 1+x. From MUSL libc (Rich Felker et al., https://musl.libc.org) src/math/log1pf.c
1858    // and fdlibm (Sun Microsystems, https://www.netlib.org/fdlibm/) s_log1p.c.
1859    // When 1+x rounds to 1 exactly (subnormal x), return x directly.
1860    // When 1+x overflows (+inf), return ln(u) without correction.
1861    // Mathematically exact: compensates for the rounding loss in 1+x without
1862    // needing a series threshold.
1863    let u = self + Self::ONE;
1864    let eq = u.simd_eq(Self::ONE);
1865    let ln_u = Self::ln(u);
1866    let correction = self * (ln_u / (u - Self::ONE));
1867    let result = eq.select(self, correction);
1868    let over = u.is_inf();
1869    over.select(ln_u, result)
1870  }
1871
1872  #[inline]
1873  pub fn sinh(self) -> Self {
1874    const_f32_as_f32x4!(P0, 1.0);
1875    const_f32_as_f32x4!(P1, 1.0 / 6.0);
1876    const_f32_as_f32x4!(P2, 1.0 / 120.0);
1877    const_f32_as_f32x4!(P3, 1.0 / 5040.0);
1878    let a = self.abs();
1879    // |x| < 0.5: Taylor poly; last truncation term < 1 ULP at x=0.5 for both types
1880    let small = a.simd_lt(f32x4::from(0.5));
1881    let t = a * a;
1882    let poly = a * polynomial_3!(t, P0, P1, P2, P3);
1883    let exp_based = {
1884      let e = a.exp();
1885      (e - Self::ONE / e) * Self::HALF
1886    };
1887    let result = small.select(poly, exp_based);
1888    result.flip_signs(self)
1889  }
1890
1891  #[inline]
1892  pub fn cosh(self) -> Self {
1893    const_f32_as_f32x4!(P0, 1.0);
1894    const_f32_as_f32x4!(P1, 1.0 / 2.0);
1895    const_f32_as_f32x4!(P2, 1.0 / 24.0);
1896    const_f32_as_f32x4!(P3, 1.0 / 720.0);
1897    let a = self.abs();
1898    // |x| < 0.5: Taylor poly; last truncation term < 1 ULP at x=0.5 for both types
1899    let small = a.simd_lt(f32x4::from(0.5));
1900    let t = a * a;
1901    let poly = polynomial_3!(t, P0, P1, P2, P3);
1902    let exp_based = {
1903      let e = a.exp();
1904      (e + Self::ONE / e) * Self::HALF
1905    };
1906    small.select(poly, exp_based)
1907  }
1908
1909  #[inline]
1910  pub fn tanh(self) -> Self {
1911    // |x| < 2e-4: tanh(x) ≈ x, error x³/3 < 16·ULP(x)
1912    // bound: x² < 48·2⁻²³ → x < 2.39e-3; 2e-4 has 10× margin
1913    // |x| > 9.011: tanh(x) = ±1 to f32 precision (e⁻²ˣ < 2⁻²⁴)
1914    let a = self.abs();
1915    let large = a.simd_gt(f32x4::from(9.011));
1916    if large.all() {
1917      return Self::ONE.flip_signs(self);
1918    }
1919    let small = a.simd_lt(f32x4::from(2e-4));
1920    let exp_based = {
1921      let t = (Self::from(-2.0) * a).exp_m1();
1922      let pos = -t / (t + Self::from(2.0));
1923      pos.flip_signs(self)
1924    };
1925    let result = small.select(self, exp_based);
1926    large.select(Self::ONE.flip_signs(self), result)
1927  }
1928}
1929
1930/// The following functionality exists only for [`f32x4`], or only for
1931/// particular types inconsistently.
1932impl f32x4 {
1933  #[inline]
1934  fn vm_pow2n(self) -> Self {
1935    const_f32_as_f32x4!(pow2_23, 8388608.0);
1936    const_f32_as_f32x4!(bias, 127.0);
1937    let a = self + (bias + pow2_23);
1938    let c = cast::<_, i32x4>(a) << 23;
1939    let std_result = cast::<_, f32x4>(c);
1940
1941    let min_exp = f32x4::from(-126.0);
1942    let is_sub = self.simd_lt(min_exp);
1943    if is_sub.any() {
1944      let valid = self.simd_ge(f32x4::from(-149.0));
1945      let shift_f = self + f32x4::from(149.0);
1946      let mut shift_i = shift_f.trunc_int();
1947      shift_i = cast::<_, i32x4>(valid).select(shift_i, i32x4::ZERO);
1948      let mantissa = i32x4::ONE << shift_i;
1949      let sub_result = cast::<_, f32x4>(mantissa);
1950      let sub_result = valid.select(sub_result, f32x4::ZERO);
1951      is_sub.select(sub_result, std_result)
1952    } else {
1953      std_result
1954    }
1955  }
1956
1957  #[inline]
1958  fn exponent(self) -> f32x4 {
1959    const_f32_as_f32x4!(pow2_23, 8388608.0);
1960    const_f32_as_f32x4!(bias, 127.0);
1961    let a = cast::<_, u32x4>(self);
1962    let b = a >> 23;
1963    let c = b | cast::<_, u32x4>(pow2_23);
1964    let d = cast::<_, f32x4>(c);
1965    let e = d - (pow2_23 + bias);
1966    e
1967  }
1968
1969  #[inline]
1970  fn fraction_2(self) -> Self {
1971    let t1 = cast::<_, u32x4>(self);
1972    let t2 = cast::<_, u32x4>(
1973      (t1 & u32x4::from(0x007FFFFF)) | u32x4::from(0x3F000000),
1974    );
1975    cast::<_, f32x4>(t2)
1976  }
1977  #[inline]
1978  fn is_zero_or_subnormal(self) -> Self {
1979    let t = cast::<_, i32x4>(self);
1980    let t = t & i32x4::splat(0x7F800000);
1981    let mask = t.simd_eq(i32x4::splat(0));
1982    cast::<_, f32x4>(mask)
1983  }
1984  #[inline]
1985  fn infinity() -> Self {
1986    cast::<_, f32x4>(i32x4::splat(0x7F800000))
1987  }
1988  #[inline]
1989  fn nan_log() -> Self {
1990    cast::<_, f32x4>(i32x4::splat(0x7FC00000 | 0x101 & 0x003FFFFF))
1991  }
1992  #[inline]
1993  fn nan_pow() -> Self {
1994    cast::<_, f32x4>(i32x4::splat(0x7FC00000 | 0x101 & 0x003FFFFF))
1995  }
1996
1997  /// Returns `[self[0], b[0], self[1], b[1]]`.
1998  #[must_use]
1999  #[inline]
2000  pub fn unpack_lo(self, b: Self) -> Self {
2001    pick! {
2002      if #[cfg(target_feature="sse")] {
2003        Self { sse: unpack_low_m128(self.sse, b.sse) }
2004      } else if #[cfg(target_feature="simd128")] {
2005        Self {
2006          simd: u32x4_shuffle::<0, 4, 1, 5>(self.simd, b.simd)
2007        }
2008      } else if #[cfg(all(target_feature="neon", target_arch="aarch64"))]{
2009        unsafe {Self { neon: vzip1q_f32(self.neon, b.neon) }}
2010      } else {
2011        Self { arr: [
2012          self.arr[0],
2013          b.arr[0],
2014          self.arr[1],
2015          b.arr[1],
2016        ]}
2017      }
2018    }
2019  }
2020
2021  /// Returns `[self[2], b[2], self[3], b[3]]`.
2022  #[must_use]
2023  #[inline]
2024  pub fn unpack_hi(self, b: Self) -> Self {
2025    pick! {
2026      if #[cfg(target_feature="sse")] {
2027        Self { sse: unpack_high_m128(self.sse, b.sse) }
2028      } else if #[cfg(target_feature="simd128")] {
2029        Self {
2030          simd: u32x4_shuffle::<2, 6, 3, 7>(self.simd, b.simd)
2031        }
2032      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))]{
2033        unsafe {Self { neon: vzip2q_f32(self.neon, b.neon) }}
2034      } else {
2035        Self { arr: [
2036          self.arr[2],
2037          b.arr[2],
2038          self.arr[3],
2039          b.arr[3],
2040        ]}
2041      }
2042    }
2043  }
2044
2045  /// Converts each element from [`i32`] to [`f32`].
2046  #[inline]
2047  pub fn from_i32x4(v: i32x4) -> Self {
2048    pick! {
2049      if #[cfg(target_feature="sse2")] {
2050        Self { sse: convert_to_m128_from_i32_m128i(v.sse) }
2051      } else if #[cfg(target_feature="simd128")] {
2052        Self { simd: f32x4_convert_i32x4(v.simd) }
2053      } else if #[cfg(all(target_feature="neon",target_arch="aarch64"))] {
2054        Self { neon: unsafe { vcvtq_f32_s32(v.neon) }}
2055      } else {
2056        Self { arr: [
2057            v.as_array()[0] as f32,
2058            v.as_array()[1] as f32,
2059            v.as_array()[2] as f32,
2060            v.as_array()[3] as f32,
2061          ] }
2062      }
2063    }
2064  }
2065
2066  /// Returns a [mask] that checks if each element has a negative sign,
2067  /// including `-0.0`, NaNs with negative sign bit and negative infinity.
2068  ///
2069  /// Note that this function has a misleading name. If the sign bit is set, the
2070  /// result has all bits set, not just the sign bit. This function has been
2071  /// renamed to [`is_sign_negative`].
2072  ///
2073  /// [mask]: crate#masks
2074  /// [`is_sign_negative`]: Self::is_sign_negative
2075  #[inline]
2076  #[must_use]
2077  #[deprecated(since = "1.4.0", note = "renamed to `is_sign_negative`")]
2078  pub fn sign_bit(self) -> Self {
2079    self.is_sign_negative()
2080  }
2081}
wide/f32x4_.rs

wide/
f32x4_.rs
