1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
//! `mul_f64x4` — packed `f64x4` multiplication.
//!
//! * `x86_64` + AVX2: `vmulpd ymm` (256-bit, 4 lanes).
//! * `AArch64`: two `fmul v.2d` NEON ops (NEON mandatory).
//! * riscv64 + RVV: `vfmul.vv` with `vsetvli` for 4×f64.
//! * fallback: scalar loop.
#![allow(unsafe_code)]
/// Multiplies two 4-lane `f64` vectors element-wise.
///
/// # Safety
///
/// On the AVX2 path uses raw pointers + 256-bit unaligned ops. Lengths
/// are checked; alignment isn't (we emit `vmovupd`).
#[allow(clippy::inline_always)]
#[inline(always)]
pub fn apply(lhs: &[f64], rhs: &[f64], out: &mut [f64]) {
debug_assert!(
lhs.len() == 4 && rhs.len() == 4 && out.len() == 4,
"mul_f64x4::apply requires exactly 4-element slices \
(got lhs={}, rhs={}, out={})",
lhs.len(),
rhs.len(),
out.len()
);
if lhs.len() != 4 || rhs.len() != 4 || out.len() != 4 {
return;
}
#[cfg(target_arch = "x86_64")]
{
if std::is_x86_feature_detected!("avx2") {
// SAFETY: lengths checked; AVX2 detected.
unsafe {
use core::arch::asm;
asm!(
"vmovupd ymm0, ymmword ptr [{lhs}]",
"vmovupd ymm1, ymmword ptr [{rhs}]",
"vmulpd ymm0, ymm0, ymm1",
"vmovupd ymmword ptr [{out}], ymm0",
lhs = in(reg) lhs.as_ptr(),
rhs = in(reg) rhs.as_ptr(),
out = in(reg) out.as_mut_ptr(),
out("ymm0") _,
out("ymm1") _,
options(nostack, preserves_flags),
);
}
return;
}
}
#[cfg(target_arch = "aarch64")]
{
// SAFETY: lengths checked above; NEON is mandatory on AArch64.
unsafe {
use core::arch::asm;
asm!(
"ld1 {{v0.2d}}, [{lhs}], #16",
"ld1 {{v1.2d}}, [{rhs}], #16",
"ld1 {{v2.2d}}, [{lhs}]",
"ld1 {{v3.2d}}, [{rhs}]",
"fmul v0.2d, v0.2d, v1.2d",
"fmul v2.2d, v2.2d, v3.2d",
"st1 {{v0.2d}}, [{out}], #16",
"st1 {{v2.2d}}, [{out}]",
lhs = inout(reg) lhs.as_ptr() => _,
rhs = inout(reg) rhs.as_ptr() => _,
out = inout(reg) out.as_mut_ptr() => _,
out("v0") _, out("v1") _,
out("v2") _, out("v3") _,
options(nostack, preserves_flags),
);
}
return;
}
#[cfg(all(target_arch = "riscv64", target_feature = "v"))]
{
// SAFETY: lengths checked above; RVV activated via target_feature = "v".
unsafe {
use core::arch::asm;
asm!(
"li t0, 4",
"vsetvli t0, t0, e64, m1, ta, ma",
"vle64.v v0, ({lhs})",
"vle64.v v1, ({rhs})",
"vfmul.vv v0, v0, v1",
"vse64.v v0, ({out})",
lhs = in(reg) lhs.as_ptr(),
rhs = in(reg) rhs.as_ptr(),
out = in(reg) out.as_mut_ptr(),
out("t0") _,
out("v0") _,
out("v1") _,
options(nostack),
);
}
return;
}
// Scalar fallback.
for i in 0..4 {
out[i] = lhs[i] * rhs[i];
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn vector_mul_matches_scalar() {
let a = [1.0_f64, 2.0, 3.0, 4.0];
let b = [10.0_f64, 20.0, 30.0, 40.0];
let mut out = [0.0_f64; 4];
apply(&a, &b, &mut out);
assert_eq!(out, [10.0, 40.0, 90.0, 160.0]);
}
#[test]
fn multiplying_by_zero_zeros_out() {
let a = [1.0_f64, 2.0, 3.0, 4.0];
let b = [0.0_f64; 4];
let mut out = [9.0_f64; 4];
apply(&a, &b, &mut out);
assert_eq!(out, [0.0_f64; 4]);
}
}