NEON backend for VectorF32 on AArch64: one float32x4_t per vector
(Len * Packing <= 4), with a vectorised sin/cos port of the SSE helpers.

VectorF16 gets a native float16x8_t path when the target has ARMv8.2-FP16
arithmetic. Without it VectorF16 aliases VectorF32, which converts
fp16 <-> fp32 on Load/Store.

The scalar fallback gains the batched Dot that MatrixRowMajor relies on.
Test tiers added for armv8-a and armv8.2-a+fp16.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
This commit is contained in:
Jorijn van der Graaf 2026-10-06 17:13:54 +02:00
commit 0c55e76873
Signed by: jorijnvdgraaf
GPG key ID: 2937E59CDCC1BCFB
6 changed files with 825 additions and 17 deletions

View file

@ -5,15 +5,19 @@ module;
#ifdef __x86_64
#include <immintrin.h>
#endif
#if defined(__aarch64__) && defined(__ARM_NEON)
#include <arm_neon.h>
#endif
export module Crafter.Math:VectorF16;
import std;
import :Common;
// A real VectorF16 struct is provided only when _Float16 is available AND we
// are either on an x86_64 build with AVX512FP16 (intrinsic impl below) or on a
// non-x86_64 target (scalar fallback further below). On x86_64 without
// AVX512FP16 or anywhere without _Float16, Crafter.Math:Basic aliases
// VectorF16 to VectorF32 instead.
// are on an x86_64 build with AVX512FP16 or an AArch64 NEON build with
// ARMv8.2-FP16 (intrinsic impls below), or on a target with neither x86_64 nor
// NEON (scalar fallback further below). On x86_64 without AVX512FP16, on
// AArch64 NEON without FP16 or anywhere without _Float16, Crafter.Math:Basic
// aliases VectorF16 to VectorF32 instead.
#ifdef __AVX512FP16__
namespace Crafter {
export template <std::uint8_t Len, std::uint8_t Packing>
@ -1049,7 +1053,379 @@ namespace Crafter {
}
#elif !defined(__x86_64) && defined(__FLT16_MAX__)
#elif defined(__aarch64__) && defined(__ARM_NEON) && defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC)
// AArch64 NEON with ARMv8.2-FP16 half-precision arithmetic. Built the same
// way as the NEON VectorF32: every vector is one 128-bit register, here a
// float16x8_t with eight lanes, and shuffles and blends are compile-time lane
// maps handed to __builtin_shufflevector.
namespace Crafter {
export template <std::uint8_t Len, std::uint8_t Packing>
struct VectorF16 : public VectorBase<Len, Packing, _Float16> {
template <std::uint8_t Len2, std::uint8_t Packing2>
friend struct VectorF16;
using Base = VectorBase<Len, Packing, _Float16>;
static constexpr std::uint8_t NElems = Base::AlignmentElement;
static_assert(NElems == 8, "NEON VectorF16 assumes 8-lane vectors");
private:
// Source lane for each of the eight output lanes: 0..7 pick from the
// first operand of Permute, 8..15 from the second. Every map names a
// real source lane for all eight outputs, padding included, so no
// lane of a result is ever undefined.
using LaneMap = std::array<int, 8>;
static constexpr LaneMap InPlace{0, 1, 2, 3, 4, 5, 6, 7};
template <LaneMap Lanes, typename Register>
static constexpr Register Permute(Register a, Register b) {
return __builtin_shufflevector(a, b, Lanes[0], Lanes[1], Lanes[2], Lanes[3], Lanes[4], Lanes[5], Lanes[6], Lanes[7]);
}
// Per-pack shuffle; padding lanes stay where they are.
template <std::array<std::uint8_t, Len> ShuffleValues>
static consteval LaneMap GetShuffleLanes() {
LaneMap lanes = InPlace;
for (std::uint8_t p = 0; p < Packing; ++p)
for (std::uint8_t i = 0; i < Len; ++i)
lanes[p * Len + i] = p * Len + ShuffleValues[i];
return lanes;
}
// Per-pack blend: the second operand where the mask is set, the first
// elsewhere.
template <std::array<bool, Len> BlendValues>
static consteval LaneMap GetBlendLanes() {
LaneMap lanes = InPlace;
for (std::uint8_t p = 0; p < Packing; ++p)
for (std::uint8_t i = 0; i < Len; ++i)
if (BlendValues[i]) lanes[p * Len + i] += 8;
return lanes;
}
// Moves every sub-vector to its slot in a (ToLen, ToPacking) layout,
// carrying over the components both layouts have. The remaining lanes
// are left as they were.
template <std::uint8_t ToLen, std::uint8_t ToPacking>
static consteval LaneMap GetResizeLanes() {
LaneMap lanes = InPlace;
for (std::uint8_t p = 0; p < std::min(Packing, ToPacking); ++p)
for (std::uint8_t i = 0; i < std::min(Len, ToLen); ++i)
lanes[p * ToLen + i] = p * Len + i;
return lanes;
}
// Normalize keeps one scale per sub-vector in a VectorF16<1, N*Packing>;
// this spreads the scales of input `Index` over the Len lanes of each
// of its sub-vectors. Padding lanes repeat the last scale.
template <std::uint8_t Index>
static consteval LaneMap GetSpreadLanes() {
LaneMap lanes{};
for (std::uint8_t i = 0; i < 8; ++i)
lanes[i] = Index * Packing + std::min(i / Len, Packing - 1);
return lanes;
}
// QuanternionFromEuler factor row out of the sin (first operand) and
// cos (second operand) of a VectorF16<3, Packing>: in each quaternion,
// lane `Axis` takes the sin of that axis and the other three lanes
// its cos, or the other way around when `Swap` is set.
template <std::uint8_t Axis, bool Swap>
static consteval LaneMap GetEulerLanes() {
LaneMap lanes = InPlace;
for (std::uint8_t p = 0; p < Packing; ++p)
for (std::uint8_t i = 0; i < 4; ++i)
lanes[p * 4 + i] = ((i == Axis) != Swap ? 0 : 8) + p * 3 + Axis;
return lanes;
}
// The batched reductions sum sub-vectors with faddp, which adds the
// adjacent lane pairs of two registers at once. That takes every
// sub-vector in a power-of-two slot of lanes: `Slot` lanes wide,
// `Slots` of them per register.
static constexpr std::uint8_t Slot = std::bit_ceil(Len);
static constexpr std::uint8_t Slots = 8 / Slot;
// Spreads the sub-vectors out to one slot each and zero-fills the rest
// (lane 8 is a zero from Permute's second operand). Power-of-two Len
// sits in its slots already, and whatever shares the register with it
// ends up in sums of its own.
static consteval LaneMap GetSlotLanes() {
LaneMap lanes{};
for (std::uint8_t i = 0; i < 8; ++i)
lanes[i] = (i / Slot < Packing && i % Slot < Len) ? (i / Slot) * Len + i % Slot : 8;
return lanes;
}
static constexpr float16x8_t ToSlots(float16x8_t m) {
if constexpr (Slot != Len) {
return Permute<GetSlotLanes()>(m, vdupq_n_f16(0));
} else {
return m;
}
}
// Once the faddp levels are done, a register holds the sums of `Slot`
// inputs, each in the leading Packing lanes of a group of Slots. This
// is the lane the i-th of those sums sits in.
static consteval int SumLane(std::uint8_t i) {
return (i / Packing) * Slots + i % Packing;
}
static consteval LaneMap GetSumLanes() {
LaneMap lanes = InPlace;
for (std::uint8_t i = 0; i < Slot * Packing; ++i) lanes[i] = SumLane(i);
return lanes;
}
// Puts the first `Width` values of a, then the first `Width` of b, at
// the front. `Sums` says the values still sit where the faddp levels
// left them (see SumLane) instead of already being at the front.
template <std::uint8_t Width, bool Sums>
static consteval LaneMap GetJoinLanes() {
LaneMap lanes = InPlace;
for (std::uint8_t i = 0; i < Width; ++i) {
lanes[i] = Sums ? SumLane(i) : i;
lanes[Width + i] = 8 + lanes[i];
}
return lanes;
}
// Concatenates the sums of the first `Count` registers, `Width` of
// them per register, halving the register count per step.
template <std::uint8_t Count, std::uint8_t Width, bool Sums>
static constexpr float16x8_t Join(std::array<float16x8_t, Base::BatchSize> r) {
if constexpr (Count > 1) {
for (std::uint8_t i = 0; i < Count / 2; ++i)
r[i] = Permute<GetJoinLanes<Width, Sums>()>(r[2 * i], r[2 * i + 1]);
return Join<Count / 2, 2 * Width, false>(r);
} else if constexpr (Sums) {
return Permute<GetSumLanes()>(r[0], r[0]);
} else {
return r[0];
}
}
// Collapses every sub-vector of the BatchSize registers into the sum
// of its components: result lane i*Packing+p is the sum of lanes
// [p*Len, p*Len+Len) of m[i]. Each faddp level halves the slot width
// and packs two registers into one, so log2(Slot) levels leave plain
// sums that Join then lines up.
static constexpr float16x8_t SumPacks(std::array<float16x8_t, Base::BatchSize> m) {
constexpr std::uint8_t N = Base::BatchSize;
static_assert(N % Slot == 0);
for (float16x8_t& product : m) product = ToSlots(product);
std::uint8_t count = N;
for (std::uint8_t width = Slot; width > 1; width /= 2) {
count /= 2;
for (std::uint8_t i = 0; i < count; ++i) m[i] = vpaddq_f16(m[2 * i], m[2 * i + 1]);
}
return Join<N / Slot, Slot * Packing, true>(m);
}
public:
constexpr VectorF16() = default;
constexpr VectorF16(float16x8_t vv) { this->v = vv; }
constexpr VectorF16(const _Float16* vB) { Load(vB); }
constexpr VectorF16(_Float16 val) { this->v = vdupq_n_f16(val); }
// Through memcpy because arm_neon.h types its fp16 pointers as __fp16,
// which _Float16 may not alias.
constexpr void Load(const _Float16* vB) { std::memcpy(&this->v, vB, sizeof(this->v)); }
constexpr void Store(_Float16* vB) const { std::memcpy(vB, &this->v, sizeof(this->v)); }
template<typename T>
constexpr std::array<_Float16, NElems> Store() const {
std::array<_Float16, NElems> r;
Store(r.data());
return r;
}
template <std::uint8_t BLen, std::uint8_t BPacking>
constexpr operator VectorF16<BLen, BPacking>() const {
return VectorF16<BLen, BPacking>(Permute<GetResizeLanes<BLen, BPacking>()>(this->v, this->v));
}
constexpr VectorF16<Len, Packing> operator+(VectorF16<Len, Packing> b) const { return VectorF16<Len, Packing>(vaddq_f16(this->v, b.v)); }
constexpr VectorF16<Len, Packing> operator-(VectorF16<Len, Packing> b) const { return VectorF16<Len, Packing>(vsubq_f16(this->v, b.v)); }
constexpr VectorF16<Len, Packing> operator*(VectorF16<Len, Packing> b) const { return VectorF16<Len, Packing>(vmulq_f16(this->v, b.v)); }
constexpr VectorF16<Len, Packing> operator/(VectorF16<Len, Packing> b) const { return VectorF16<Len, Packing>(vdivq_f16(this->v, b.v)); }
constexpr void operator+=(VectorF16<Len, Packing> b) { this->v = vaddq_f16(this->v, b.v); }
constexpr void operator-=(VectorF16<Len, Packing> b) { this->v = vsubq_f16(this->v, b.v); }
constexpr void operator*=(VectorF16<Len, Packing> b) { this->v = vmulq_f16(this->v, b.v); }
constexpr void operator/=(VectorF16<Len, Packing> b) { this->v = vdivq_f16(this->v, b.v); }
constexpr VectorF16<Len, Packing> operator+(_Float16 b) const { return *this + VectorF16<Len, Packing>(b); }
constexpr VectorF16<Len, Packing> operator-(_Float16 b) const { return *this - VectorF16<Len, Packing>(b); }
constexpr VectorF16<Len, Packing> operator*(_Float16 b) const { return *this * VectorF16<Len, Packing>(b); }
constexpr VectorF16<Len, Packing> operator/(_Float16 b) const { return *this / VectorF16<Len, Packing>(b); }
constexpr void operator+=(_Float16 b) { *this += VectorF16<Len, Packing>(b); }
constexpr void operator-=(_Float16 b) { *this -= VectorF16<Len, Packing>(b); }
constexpr void operator*=(_Float16 b) { *this *= VectorF16<Len, Packing>(b); }
constexpr void operator/=(_Float16 b) { *this /= VectorF16<Len, Packing>(b); }
constexpr VectorF16<Len, Packing> operator-() const { return VectorF16<Len, Packing>(vnegq_f16(this->v)); }
// Only the Len*Packing lanes in use take part; padding lanes can hold
// anything. Narrowing the compare result to a byte per lane puts all
// eight in one 64-bit integer, where the lanes in use are its low
// bytes.
constexpr bool operator==(VectorF16<Len, Packing> b) const {
constexpr std::uint64_t used = ~0ull >> (64 - 8 * Len * Packing);
std::uint64_t equal = vget_lane_u64(vreinterpret_u64_u8(vmovn_u16(vceqq_f16(this->v, b.v))), 0);
return (equal & used) == used;
}
constexpr bool operator!=(VectorF16<Len, Packing> b) const { return !(*this == b); }
template<std::uint32_t ExtractLen>
constexpr VectorF16<ExtractLen, Packing> ExtractLo() const {
return VectorF16<ExtractLen, Packing>(Permute<GetResizeLanes<ExtractLen, Packing>()>(this->v, this->v));
}
// Sin and cos run in single precision: the register widens to fp32,
// goes through the shared kernels and narrows back. The upper half is
// skipped when it is all padding.
constexpr std::tuple<VectorF16<Len, Packing>, VectorF16<Len, Packing>> SinCos() const {
float32x4_t sinLo, cosLo;
Base::sincos_f32x4(vcvt_f32_f16(vget_low_f16(this->v)), sinLo, cosLo);
float32x4_t sinHi = sinLo, cosHi = cosLo;
if constexpr (Len * Packing > 4) {
Base::sincos_f32x4(vcvt_f32_f16(vget_high_f16(this->v)), sinHi, cosHi);
}
return {
VectorF16<Len, Packing>(vcombine_f16(vcvt_f16_f32(sinLo), vcvt_f16_f32(sinHi))),
VectorF16<Len, Packing>(vcombine_f16(vcvt_f16_f32(cosLo), vcvt_f16_f32(cosHi)))
};
}
constexpr VectorF16<Len, Packing> Cos() const { return std::get<1>(SinCos()); }
constexpr VectorF16<Len, Packing> Sin() const { return std::get<0>(SinCos()); }
template <std::array<bool, Len> values>
constexpr VectorF16<Len, Packing> Negate() const {
constexpr auto mask = []() {
std::array<std::uint16_t, NElems> m{};
for (std::uint8_t p = 0; p < Packing; ++p)
for (std::uint8_t i = 0; i < Len; ++i)
m[p * Len + i] = values[i] ? 0x8000u : 0u;
return m;
}();
return VectorF16<Len, Packing>(vreinterpretq_f16_u16(veorq_u16(vreinterpretq_u16_f16(this->v), vld1q_u16(mask.data()))));
}
// vfmaq_f16(a, b, c) is a + b*c.
static constexpr VectorF16<Len, Packing> MulitplyAdd(VectorF16<Len, Packing> a, VectorF16<Len, Packing> b, VectorF16<Len, Packing> add) {
return VectorF16<Len, Packing>(vfmaq_f16(add.v, a.v, b.v));
}
static constexpr VectorF16<Len, Packing> MulitplySub(VectorF16<Len, Packing> a, VectorF16<Len, Packing> b, VectorF16<Len, Packing> sub) {
return VectorF16<Len, Packing>(vfmaq_f16(vnegq_f16(sub.v), a.v, b.v));
}
// a × b = (a * b.yzx - b * a.yzx).yzx, one shuffle less than building
// both yzx and zxy of each input.
constexpr static VectorF16<Len, Packing> Cross(VectorF16<Len, Packing> a, VectorF16<Len, Packing> b) requires(Len == 3) {
VectorF16<Len, Packing> aYzx = a.template Shuffle<{{1,2,0}}>();
VectorF16<Len, Packing> bYzx = b.template Shuffle<{{1,2,0}}>();
VectorF16<Len, Packing> zxy(vfmaq_f16(vmulq_f16(a.v, bYzx.v), vnegq_f16(b.v), aYzx.v));
return zxy.template Shuffle<{{1,2,0}}>();
}
constexpr static _Float16 Dot(VectorF16<Len, Packing> a, VectorF16<Len, Packing> b) requires(Packing == 1) {
float16x8_t sum = ToSlots(vmulq_f16(a.v, b.v));
for (std::uint8_t width = Slot; width > 1; width /= 2) sum = vpaddq_f16(sum, sum);
return vgetq_lane_f16(sum, 0);
}
template <const std::array<std::uint8_t, Len> ShuffleValues>
constexpr VectorF16<Len, Packing> Shuffle() const {
return VectorF16<Len, Packing>(Permute<GetShuffleLanes<ShuffleValues>()>(this->v, this->v));
}
template <std::array<bool, Len> ShuffleValues>
constexpr static VectorF16<Len, Packing> Blend(VectorF16<Len, Packing> a, VectorF16<Len, Packing> b) {
return VectorF16<Len, Packing>(Permute<GetBlendLanes<ShuffleValues>()>(a.v, b.v));
}
template<typename... Rest>
requires((std::is_same_v<Rest, VectorF16<Len, Packing>> && ...) &&
(1 + sizeof...(Rest) == VectorBase<Len, Packing, _Float16>::BatchSize))
constexpr static auto LengthSq(VectorF16<Len, Packing> first, Rest... rest) {
constexpr std::uint8_t N = VectorBase<Len, Packing, _Float16>::BatchSize;
std::array<VectorF16<Len, Packing>, N> args{ first, rest... };
std::array<float16x8_t, N> squares;
for (std::uint8_t i = 0; i < N; ++i) squares[i] = vmulq_f16(args[i].v, args[i].v);
return VectorF16<1, static_cast<std::uint8_t>(Packing * N)>(SumPacks(squares));
}
template<typename... Rest>
requires((std::is_same_v<Rest, VectorF16<Len, Packing>> && ...) &&
(1 + sizeof...(Rest) == VectorBase<Len, Packing, _Float16>::BatchSize))
constexpr static auto Length(VectorF16<Len, Packing> first, Rest... rest) {
auto sq = LengthSq(first, rest...);
sq.v = vsqrtq_f16(sq.v);
return sq;
}
// Pairwise dot products across BatchSize pairs: Dot(a0, b0, a1, b1,
// ...) returns one lane per sub-vector pair.
template<typename... Rest>
requires((std::is_same_v<Rest, VectorF16<Len, Packing>> && ...) &&
(1 + sizeof...(Rest) == 2 * VectorBase<Len, Packing, _Float16>::BatchSize))
constexpr static auto Dot(VectorF16<Len, Packing> first, Rest... rest) {
constexpr std::uint8_t N = VectorBase<Len, Packing, _Float16>::BatchSize;
std::array<VectorF16<Len, Packing>, 2 * N> args{ first, rest... };
std::array<float16x8_t, N> products;
for (std::uint8_t i = 0; i < N; ++i) products[i] = vmulq_f16(args[2 * i].v, args[2 * i + 1].v);
return VectorF16<1, static_cast<std::uint8_t>(Packing * N)>(SumPacks(products));
}
// One sqrt and one divide cover the whole batch; each input then
// multiplies by its own lanes of the result.
template<typename... Rest>
requires((std::is_same_v<Rest, VectorF16<Len, Packing>> && ...) &&
(1 + sizeof...(Rest) == VectorBase<Len, Packing, _Float16>::BatchSize))
constexpr static auto Normalize(VectorF16<Len, Packing> first, Rest... rest) {
constexpr std::uint8_t N = VectorBase<Len, Packing, _Float16>::BatchSize;
std::array<VectorF16<Len, Packing>, N> args{ first, rest... };
float16x8_t scale = vdivq_f16(vdupq_n_f16(1), Length(first, rest...).v);
return [&]<std::size_t... Is>(std::index_sequence<Is...>) {
return std::array<VectorF16<Len, Packing>, N>{
VectorF16<Len, Packing>(vmulq_f16(args[Is].v, Permute<GetSpreadLanes<Is>()>(scale, scale)))...
};
}(std::make_index_sequence<N>{});
}
constexpr static VectorF16<Len, Packing> Rotate(VectorF16<3, Packing> v, VectorF16<4, Packing> q) requires(Len == 3) {
VectorF16<3, Packing> qv = q;
VectorF16<3, Packing> qw = q.template Shuffle<{{3,3,3,3}}>();
VectorF16<3, Packing> t = Cross(qv, v) * _Float16(2);
return MulitplyAdd(t, qw, v) + Cross(qv, t);
}
constexpr static VectorF16<3, Packing> RotatePivot(VectorF16<3, Packing> v, VectorF16<4, Packing> q, VectorF16<3, Packing> pivot) requires(Len == 3) {
VectorF16<3, Packing> translated = v - pivot;
return Rotate(translated, q) + pivot;
}
// With (roll, pitch, yaw) half angles in x, y, z:
// x = sx*cy*cz - cx*sy*sz y = cx*sy*cz + sx*cy*sz
// z = cx*cy*sz - sx*sy*cz w = cx*cy*cz + sx*sy*sz
// Each factor row below is one two-source shuffle of the sin and cos
// registers (see GetEulerLanes).
constexpr static VectorF16<4, Packing> QuanternionFromEuler(VectorF16<3, Packing> EulerHalf) requires(Len == 4) {
std::tuple<VectorF16<3, Packing>, VectorF16<3, Packing>> sinCos = EulerHalf.SinCos();
float16x8_t s = std::get<0>(sinCos).v;
float16x8_t c = std::get<1>(sinCos).v;
float16x8_t left = vmulq_f16(
vmulq_f16(
Permute<GetEulerLanes<0, false>()>(s, c), // sx cx cx cx
Permute<GetEulerLanes<1, false>()>(s, c)), // cy sy cy cy
Permute<GetEulerLanes<2, false>()>(s, c)); // cz cz sz cz
float16x8_t right = vmulq_f16(
Permute<GetEulerLanes<0, true>()>(s, c), // cx sx sx sx
Permute<GetEulerLanes<1, true>()>(s, c)); // sy cy sy sy
VectorF16<4, Packing> rightZ(Permute<GetEulerLanes<2, true>()>(s, c)); // sz sz cz sz
rightZ = rightZ.template Negate<{{true, false, true, false}}>();
return VectorF16<4, Packing>(vfmaq_f16(left, right, rightZ.v));
}
};
}
#elif !defined(__x86_64) && defined(__FLT16_MAX__) && !(defined(__aarch64__) && defined(__ARM_NEON))
// Scalar software fallback for non-x86_64 targets that still have _Float16.
namespace Crafter {
export template <std::uint8_t Len, std::uint8_t Packing>
@ -1304,7 +1680,7 @@ namespace Crafter {
}
#endif
#if defined(__FLT16_MAX__) && (!defined(__x86_64) || defined(__AVX512FP16__))
#if defined(__FLT16_MAX__) && (!defined(__x86_64) || defined(__AVX512FP16__)) && !(defined(__aarch64__) && defined(__ARM_NEON) && !defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC))
export template <std::uint32_t Len, std::uint32_t Packing>
struct std::formatter<Crafter::VectorF16<Len, Packing>> : std::formatter<std::string> {
constexpr auto format(const Crafter::VectorF16<Len, Packing>& obj, format_context& ctx) const {