From 0c55e76873ae4cb08f3544d30267ede89462d0ea Mon Sep 17 00:00:00 2001 From: Jorijn van der Graaf Date: Tue, 6 Oct 2026 17:13:54 +0200 Subject: [PATCH] ARM NEON 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 --- README.md | 9 +- interfaces/Crafter.Math-Basic.cppm | 2 +- interfaces/Crafter.Math-Common.cppm | 94 +++++- interfaces/Crafter.Math-VectorF16.cppm | 388 ++++++++++++++++++++++++- interfaces/Crafter.Math-VectorF32.cppm | 339 ++++++++++++++++++++- project.cpp | 10 + 6 files changed, 825 insertions(+), 17 deletions(-) diff --git a/README.md b/README.md index 04f72ee..9135324 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ A C++23 math library for the Crafter engine, distributed as a set of C++ modules. Provides generic vector and matrix types alongside SIMD-specialized fixed-size vectors for `float` and `_Float16`, plus a small set of ray intersection routines. -The library is hardware-aware: it picks the widest SIMD path the target supports (SSE / AVX / AVX-512 on x86_64, RVV on RISC-V, SIMD128 on WebAssembly) and falls back to scalar code elsewhere. +The library is hardware-aware: it picks the widest SIMD path the target supports (SSE / AVX / AVX-512 on x86_64, NEON on AArch64, RVV on RISC-V, SIMD128 on WebAssembly) and falls back to scalar code elsewhere. ## Modules @@ -10,13 +10,13 @@ The library is hardware-aware: it picks the widest SIMD path the target supports | Partition | Contents | | --- | --- | -| `:Basic` | `ToRadian`, plus the `VectorF16 → VectorF32` alias on x86_64 without AVX512-FP16. | +| `:Basic` | `ToRadian`, plus the `VectorF16 → VectorF32` alias on x86_64 without AVX512-FP16 and on AArch64 without ARMv8.2-FP16. | | `:Vector` | Generic `Vector` with `.x/.y/.z/.w` and `.r/.g/.b/.a` accessors, arithmetic, comparison, and conversion operators. | | `:MatrixRowMajor` | `MatrixRowMajor` — row-major matrices with packed-batch support via `Repeats`. | | `:VectorF32` | `VectorF32` — SIMD-specialized 32-bit float vector with `Load`/`Store`, FMA, `Cross`, `Dot`, `Length`, `Normalize`, `Sin`/`Cos`, `Shuffle`, etc. | | `:VectorF16` | `VectorF16` — same surface as `VectorF32` for `_Float16` when the target supports it. | | `:Intersection` | `IntersectionTestRayTriangle`, `IntersectionTestRaySphere`, `IntersectionTestRayOrientedBox`. | -| `:Common` | Internal SIMD type selection (`__m128/256/512`, `__m128h/256h/512h`, `v128_t`, or scalar array) used by the F16/F32 vectors. | +| `:Common` | Internal SIMD type selection (`__m128/256/512`, `__m128h/256h/512h`, `float32x4_t`, `float16x8_t`, `v128_t`, or scalar array) used by the F16/F32 vectors. | `Len` is the logical vector width (1–4 for most graphics use) and `Packing` is the number of vectors packed into one SIMD register — e.g. `VectorF32<3, 4>` holds four 3-vectors in a single 512-bit register, which is what the batched `Dot`, `Length`, and `Normalize` overloads operate on. @@ -52,6 +52,8 @@ This project uses [Crafter.Build](https://forgejo.catcrafts.net/Catcrafts/Crafte - `Vector-rv64gcv_zvl512b` — RVV with VLEN≥512 (uses the 64-byte tier) - `Vector-rv64gcv_zvl256b` — RVV with VLEN≥256 (32-byte tier) - `Vector-rv64gcv` — RVA23 baseline ZVL128B (16-byte tier) +- `Vector-armv8.2-a+fp16` — AArch64 NEON with ARMv8.2-FP16 (native `VectorF16`) +- `Vector-armv8-a` — AArch64 baseline NEON Build and run tests via your usual Crafter.Build entry point; `build/` and `bin/` are git-ignored. @@ -59,6 +61,7 @@ Build and run tests via your usual Crafter.Build entry point; `build/` and `bin/ - **x86_64**: SSE / AVX / AVX-512F selected per target; AVX512-FP16 is required for native `VectorF16` (otherwise it aliases `VectorF32`). F16C is used for fp16 ↔ fp32 conversion when available. - **RISC-V**: RVV 1.0 (the `V` extension, ratified Nov 2021, mandatory in the RVA23 profile). Storage is a 16/32/64-byte GNU vector picked at compile time from the guaranteed VLEN (`__riscv_v_fixed_vlen` from Clang's `-mrvv-vector-bits=`, else `__riscv_v_min_vlen` from the march's `Zvl*` suffix, else the ZVL128B baseline). `VectorF16` aliases `VectorF32` until a `Zvfh` path lands. +- **AArch64**: NEON (`arm_neon.h`) 128-bit path; SIMD width is capped at 16 bytes, so a vector holds at most 4 floats or 8 halves. `VectorF32` only needs the ARMv8-A baseline. ARMv8.2-FP16 (`-march=armv8.2-a+fp16`, or any `-march`/`-mcpu` that implies it) is required for native `VectorF16`; otherwise it aliases `VectorF32`, which converts fp16 ↔ fp32 on `Load`/`Store`. 32-bit ARM takes the scalar fallback. - **WebAssembly**: `wasm_simd128.h` 128-bit path; SIMD width is capped at 16 bytes. - **Other targets**: scalar fallback via `std::array`. diff --git a/interfaces/Crafter.Math-Basic.cppm b/interfaces/Crafter.Math-Basic.cppm index cd16528..7c3b3d0 100755 --- a/interfaces/Crafter.Math-Basic.cppm +++ b/interfaces/Crafter.Math-Basic.cppm @@ -11,7 +11,7 @@ namespace Crafter { return degrees * (std::numbers::pi / 180); } - #if (defined(__x86_64) && !defined(__AVX512FP16__)) || defined(__riscv_vector) || !defined(__FLT16_MAX__) + #if (defined(__x86_64) && !defined(__AVX512FP16__)) || defined(__riscv_vector) || (defined(__aarch64__) && defined(__ARM_NEON) && !defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC)) || !defined(__FLT16_MAX__) export template using VectorF16 = VectorF32; #endif diff --git a/interfaces/Crafter.Math-Common.cppm b/interfaces/Crafter.Math-Common.cppm index 6912e91..00add39 100644 --- a/interfaces/Crafter.Math-Common.cppm +++ b/interfaces/Crafter.Math-Common.cppm @@ -8,6 +8,9 @@ module; #ifdef __wasm_simd128__ #include #endif +#if defined(__aarch64__) && defined(__ARM_NEON) +#include +#endif #ifdef __riscv_vector #include // Compile-time VLEN selection. RVV is VLA at the ISA level, but storage in @@ -37,9 +40,11 @@ import std; // VectorF16 exists as a real struct when _Float16 is available AND we are not // on x86_64 without AVX512FP16 (that path aliases VectorF16 to VectorF32 in // Crafter.Math:Basic for performance). The same alias kicks in on RISC-V until -// a Zvfh path lands. Each translation unit that needs this distinction -// redefines the same condition since macros do not cross module boundaries. -#if defined(__FLT16_MAX__) && (!defined(__x86_64) || defined(__AVX512FP16__)) && !defined(__riscv_vector) +// a Zvfh path lands, and on AArch64 NEON without ARMv8.2-FP16 arithmetic, +// where VectorF32 converts fp16 <-> fp32 on Load/Store. Each translation unit +// that needs this distinction redefines the same condition since macros do not +// cross module boundaries. +#if defined(__FLT16_MAX__) && (!defined(__x86_64) || defined(__AVX512FP16__)) && !defined(__riscv_vector) && !(defined(__aarch64__) && defined(__ARM_NEON) && !defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC)) namespace Crafter { export template struct VectorF16; @@ -52,7 +57,7 @@ namespace Crafter { template struct VectorBase { - #if defined(__FLT16_MAX__) && (!defined(__x86_64) || defined(__AVX512FP16__)) && !defined(__riscv_vector) + #if defined(__FLT16_MAX__) && (!defined(__x86_64) || defined(__AVX512FP16__)) && !defined(__riscv_vector) && !(defined(__aarch64__) && defined(__ARM_NEON) && !defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC)) template friend struct VectorF16; #endif @@ -89,6 +94,11 @@ namespace Crafter { >; #elif defined(__wasm_simd128__) using VectorType = v128_t; + #elif defined(__aarch64__) && defined(__ARM_NEON) + // NEON registers are 128 bits wide, so every instantiation is one + // register: four floats, or eight halves for the native VectorF16 + // that exists when the target has ARMv8.2-FP16 arithmetic. + using VectorType = std::conditional_t, float32x4_t, float16x8_t>; #elif defined(__riscv_vector) // RVV tier mirrors the x86 selector: pick the widest register the // toolchain guarantees, then size each instantiation down to the @@ -118,6 +128,10 @@ namespace Crafter { // WASM SIMD only has 128-bit vectors; cap at 16 bytes so the entire // VectorType always fits in a single v128_t. static constexpr std::uint8_t Max = 16; + #elif defined(__aarch64__) && defined(__ARM_NEON) + // NEON only has 128-bit vectors; cap at 16 bytes so the entire + // VectorType always fits in a single register. + static constexpr std::uint8_t Max = 16; #elif defined(__riscv_vector) // RVV tier selected at compile time from the guaranteed VLEN. ZVL128B // is the RVA23 baseline; ZVL256B / ZVL512B unlock wider registers @@ -698,5 +712,77 @@ namespace Crafter { out_sin = _mm512_xor_ps(out_sin, x_sign); } #endif + + #if defined(__aarch64__) && defined(__ARM_NEON) + // --- 128-bit (NEON) helpers, same algorithm as the SSE ones above --- + // Reduce |x| into [-pi/4, pi/4], return reduced value and quadrant + static constexpr void range_reduce_f32x4(float32x4_t ax, float32x4_t& r, float32x4_t& r2, int32x4_t& q) { + float32x4_t fq = vrndnq_f32(vmulq_f32(ax, vdupq_n_f32(two_over_pi))); + q = vcvtq_s32_f32(fq); + // ax - fq * pi/2 as an fmla with the constant negated, not as + // vfmsq_f32: that one is a wrapper function in clang's arm_neon.h, + // and calling it from two module units (here and :VectorF32) fails + // some importers with "definition with same mangled name as + // another definition". + r = vfmaq_f32(ax, fq, vdupq_n_f32(-pi_over_2_hi)); + r = vfmaq_f32(r, fq, vdupq_n_f32(-pi_over_2_lo)); + r2 = vmulq_f32(r, r); + } + + // vfmaq_f32(a, b, c) is a + b*c, so each Horner step reads + // "next coefficient + accumulator * r2". + static constexpr void sincos_poly_f32x4(float32x4_t r, float32x4_t r2, float32x4_t& cos_r, float32x4_t& sin_r) { + cos_r = vfmaq_f32(vdupq_n_f32(c8), vdupq_n_f32(c10), r2); + cos_r = vfmaq_f32(vdupq_n_f32(c6), cos_r, r2); + cos_r = vfmaq_f32(vdupq_n_f32(c4), cos_r, r2); + cos_r = vfmaq_f32(vdupq_n_f32(c2), cos_r, r2); + cos_r = vfmaq_f32(vdupq_n_f32(c0), cos_r, r2); + + sin_r = vfmaq_f32(vdupq_n_f32(s7), vdupq_n_f32(s9), r2); + sin_r = vfmaq_f32(vdupq_n_f32(s5), sin_r, r2); + sin_r = vfmaq_f32(vdupq_n_f32(s3), sin_r, r2); + sin_r = vfmaq_f32(vdupq_n_f32(s1), sin_r, r2); + sin_r = vfmaq_f32(vdupq_n_f32(1.0f), sin_r, r2); + sin_r = vmulq_f32(sin_r, r); + } + + static constexpr void sincos_f32x4(float32x4_t x, float32x4_t& out_sin, float32x4_t& out_cos) { + uint32x4_t x_sign = vandq_u32(vreinterpretq_u32_f32(x), vdupq_n_u32(0x80000000u)); + float32x4_t ax = vabsq_f32(x); + + float32x4_t r, r2; int32x4_t q; + range_reduce_f32x4(ax, r, r2, q); + + float32x4_t cos_r, sin_r; + sincos_poly_f32x4(r, r2, cos_r, sin_r); + + uint32x4_t is_odd = vtstq_s32(q, vdupq_n_s32(1)); + + // cos: swap on odd, negate if (q+1)&2 + out_cos = vbslq_f32(is_odd, sin_r, cos_r); + int32x4_t cos_neg = vandq_s32(vaddq_s32(q, vdupq_n_s32(1)), vdupq_n_s32(2)); + out_cos = vreinterpretq_f32_s32(veorq_s32(vreinterpretq_s32_f32(out_cos), vshlq_n_s32(cos_neg, 30))); + + // sin: swap on odd, negate if q&2, apply input sign + out_sin = vbslq_f32(is_odd, cos_r, sin_r); + int32x4_t sin_neg = vandq_s32(q, vdupq_n_s32(2)); + out_sin = vreinterpretq_f32_s32(veorq_s32(vreinterpretq_s32_f32(out_sin), vshlq_n_s32(sin_neg, 30))); + out_sin = vreinterpretq_f32_u32(veorq_u32(vreinterpretq_u32_f32(out_sin), x_sign)); + } + + // The single-result variants share sincos_f32x4; the half that is not + // returned is dead code the compiler drops after inlining. + static constexpr float32x4_t cos_f32x4(float32x4_t x) { + float32x4_t s, c; + sincos_f32x4(x, s, c); + return c; + } + + static constexpr float32x4_t sin_f32x4(float32x4_t x) { + float32x4_t s, c; + sincos_f32x4(x, s, c); + return s; + } + #endif }; } \ No newline at end of file diff --git a/interfaces/Crafter.Math-VectorF16.cppm b/interfaces/Crafter.Math-VectorF16.cppm index 13f0bc9..e0c92ab 100755 --- a/interfaces/Crafter.Math-VectorF16.cppm +++ b/interfaces/Crafter.Math-VectorF16.cppm @@ -5,15 +5,19 @@ module; #ifdef __x86_64 #include #endif +#if defined(__aarch64__) && defined(__ARM_NEON) +#include +#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 @@ -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 + struct VectorF16 : public VectorBase { + template + friend struct VectorF16; + using Base = VectorBase; + 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; + static constexpr LaneMap InPlace{0, 1, 2, 3, 4, 5, 6, 7}; + + template + 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 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 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 + 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 + 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 + 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(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 + 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 + static constexpr float16x8_t Join(std::array r) { + if constexpr (Count > 1) { + for (std::uint8_t i = 0; i < Count / 2; ++i) + r[i] = Permute()>(r[2 * i], r[2 * i + 1]); + return Join(r); + } else if constexpr (Sums) { + return Permute(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 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(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 + constexpr std::array<_Float16, NElems> Store() const { + std::array<_Float16, NElems> r; + Store(r.data()); + return r; + } + + template + constexpr operator VectorF16() const { + return VectorF16(Permute()>(this->v, this->v)); + } + + constexpr VectorF16 operator+(VectorF16 b) const { return VectorF16(vaddq_f16(this->v, b.v)); } + constexpr VectorF16 operator-(VectorF16 b) const { return VectorF16(vsubq_f16(this->v, b.v)); } + constexpr VectorF16 operator*(VectorF16 b) const { return VectorF16(vmulq_f16(this->v, b.v)); } + constexpr VectorF16 operator/(VectorF16 b) const { return VectorF16(vdivq_f16(this->v, b.v)); } + constexpr void operator+=(VectorF16 b) { this->v = vaddq_f16(this->v, b.v); } + constexpr void operator-=(VectorF16 b) { this->v = vsubq_f16(this->v, b.v); } + constexpr void operator*=(VectorF16 b) { this->v = vmulq_f16(this->v, b.v); } + constexpr void operator/=(VectorF16 b) { this->v = vdivq_f16(this->v, b.v); } + + constexpr VectorF16 operator+(_Float16 b) const { return *this + VectorF16(b); } + constexpr VectorF16 operator-(_Float16 b) const { return *this - VectorF16(b); } + constexpr VectorF16 operator*(_Float16 b) const { return *this * VectorF16(b); } + constexpr VectorF16 operator/(_Float16 b) const { return *this / VectorF16(b); } + constexpr void operator+=(_Float16 b) { *this += VectorF16(b); } + constexpr void operator-=(_Float16 b) { *this -= VectorF16(b); } + constexpr void operator*=(_Float16 b) { *this *= VectorF16(b); } + constexpr void operator/=(_Float16 b) { *this /= VectorF16(b); } + + constexpr VectorF16 operator-() const { return VectorF16(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 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 b) const { return !(*this == b); } + + template + constexpr VectorF16 ExtractLo() const { + return VectorF16(Permute()>(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> 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(vcombine_f16(vcvt_f16_f32(sinLo), vcvt_f16_f32(sinHi))), + VectorF16(vcombine_f16(vcvt_f16_f32(cosLo), vcvt_f16_f32(cosHi))) + }; + } + constexpr VectorF16 Cos() const { return std::get<1>(SinCos()); } + constexpr VectorF16 Sin() const { return std::get<0>(SinCos()); } + + template values> + constexpr VectorF16 Negate() const { + constexpr auto mask = []() { + std::array 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(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 MulitplyAdd(VectorF16 a, VectorF16 b, VectorF16 add) { + return VectorF16(vfmaq_f16(add.v, a.v, b.v)); + } + static constexpr VectorF16 MulitplySub(VectorF16 a, VectorF16 b, VectorF16 sub) { + return VectorF16(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 Cross(VectorF16 a, VectorF16 b) requires(Len == 3) { + VectorF16 aYzx = a.template Shuffle<{{1,2,0}}>(); + VectorF16 bYzx = b.template Shuffle<{{1,2,0}}>(); + VectorF16 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 a, VectorF16 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 ShuffleValues> + constexpr VectorF16 Shuffle() const { + return VectorF16(Permute()>(this->v, this->v)); + } + + template ShuffleValues> + constexpr static VectorF16 Blend(VectorF16 a, VectorF16 b) { + return VectorF16(Permute()>(a.v, b.v)); + } + + template + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == VectorBase::BatchSize)) + constexpr static auto LengthSq(VectorF16 first, Rest... rest) { + constexpr std::uint8_t N = VectorBase::BatchSize; + std::array, N> args{ first, rest... }; + std::array 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(Packing * N)>(SumPacks(squares)); + } + + template + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == VectorBase::BatchSize)) + constexpr static auto Length(VectorF16 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 + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == 2 * VectorBase::BatchSize)) + constexpr static auto Dot(VectorF16 first, Rest... rest) { + constexpr std::uint8_t N = VectorBase::BatchSize; + std::array, 2 * N> args{ first, rest... }; + std::array 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(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 + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == VectorBase::BatchSize)) + constexpr static auto Normalize(VectorF16 first, Rest... rest) { + constexpr std::uint8_t N = VectorBase::BatchSize; + std::array, N> args{ first, rest... }; + float16x8_t scale = vdivq_f16(vdupq_n_f16(1), Length(first, rest...).v); + return [&](std::index_sequence) { + return std::array, N>{ + VectorF16(vmulq_f16(args[Is].v, Permute()>(scale, scale)))... + }; + }(std::make_index_sequence{}); + } + + constexpr static VectorF16 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>> 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()>(s, c), // sx cx cx cx + Permute()>(s, c)), // cy sy cy cy + Permute()>(s, c)); // cz cz sz cz + float16x8_t right = vmulq_f16( + Permute()>(s, c), // cx sx sx sx + Permute()>(s, c)); // sy cy sy sy + VectorF16<4, Packing> rightZ(Permute()>(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 @@ -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 struct std::formatter> : std::formatter { constexpr auto format(const Crafter::VectorF16& obj, format_context& ctx) const { diff --git a/interfaces/Crafter.Math-VectorF32.cppm b/interfaces/Crafter.Math-VectorF32.cppm index f54420f..7f832bc 100755 --- a/interfaces/Crafter.Math-VectorF32.cppm +++ b/interfaces/Crafter.Math-VectorF32.cppm @@ -11,6 +11,9 @@ module; #ifdef __riscv_vector #include #endif +#if defined(__aarch64__) && defined(__ARM_NEON) +#include +#endif export module Crafter.Math:VectorF32; import std; import :Common; @@ -1931,10 +1934,318 @@ namespace Crafter { return r; } }; +#elif defined(__aarch64__) && defined(__ARM_NEON) + // AArch64 NEON implementation. VectorType is always a float32x4_t: NEON + // registers are 128 bits wide, so Common.cppm caps Len*Packing at 4 + // floats, as on WASM. Shuffles and blends have compile-time lane maps and + // go through __builtin_shufflevector, which lets the compiler pick the + // cheapest of dup/ext/rev/zip/uzp/ins/tbl per map instead of always + // paying for a tbl and its index load. + export template + struct VectorF32 : public VectorBase { + template + friend struct VectorF32; + using Base = VectorBase; + static constexpr std::uint8_t NElems = Base::AlignmentElement; + static_assert(NElems == 4, "NEON VectorF32 assumes 4-lane vectors"); + + private: + // Source lane for each of the four output lanes: 0..3 pick from the + // first operand of Permute, 4..7 from the second. Every map names a + // real source lane for all four outputs, padding included, so no lane + // of a result is ever undefined. + using LaneMap = std::array; + + template + static constexpr float32x4_t Permute(float32x4_t a, float32x4_t b) { + return __builtin_shufflevector(a, b, Lanes[0], Lanes[1], Lanes[2], Lanes[3]); + } + + // Per-pack shuffle; padding lanes stay where they are. + template ShuffleValues> + static consteval LaneMap GetShuffleLanes() { + LaneMap lanes{0, 1, 2, 3}; + 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 BlendValues> + static consteval LaneMap GetBlendLanes() { + LaneMap lanes{0, 1, 2, 3}; + 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] += 4; + 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 + static consteval LaneMap GetResizeLanes() { + LaneMap lanes{0, 1, 2, 3}; + 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 VectorF32<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, so the map + // is a plain dup wherever it can be. + template + static consteval LaneMap GetSpreadLanes() { + LaneMap lanes{}; + for (std::uint8_t i = 0; i < 4; ++i) + lanes[i] = Index * Packing + std::min(i / Len, Packing - 1); + return lanes; + } + + // 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]. faddp adds the adjacent lane pairs of + // two registers in one go, so two levels of it reduce four 4-lane + // registers to four sums without a transpose. + static constexpr float32x4_t SumPacks(std::array m) { + if constexpr (Packing == 1 && Len >= 3) { + if constexpr (Len == 3) { + // Lane 3 is padding and must not reach the sum. + for (float32x4_t& product : m) product = vsetq_lane_f32(0.0f, product, 3); + } + return vpaddq_f32(vpaddq_f32(m[0], m[1]), vpaddq_f32(m[2], m[3])); + } else if constexpr (Packing == 1 && Len == 2) { + return vpaddq_f32( + vcombine_f32(vget_low_f32(m[0]), vget_low_f32(m[1])), + vcombine_f32(vget_low_f32(m[2]), vget_low_f32(m[3]))); + } else if constexpr (Packing == 1) { + return vcombine_f32( + vget_low_f32(vzip1q_f32(m[0], m[1])), + vget_low_f32(vzip1q_f32(m[2], m[3]))); + } else if constexpr (Packing == 2 && Len == 2) { + return vpaddq_f32(m[0], m[1]); + } else if constexpr (Packing == 2) { + return vcombine_f32(vget_low_f32(m[0]), vget_low_f32(m[1])); + } else { + // Len == 1 with a lone input: its lanes already are the sums. + return m[0]; + } + } + + public: + constexpr VectorF32() = default; + constexpr VectorF32(float32x4_t vv) { this->v = vv; } + constexpr VectorF32(const float* vB) { Load(vB); } + constexpr VectorF32(const _Float16* vB) { Load(vB); } + constexpr VectorF32(float val) { this->v = vdupq_n_f32(val); } + + constexpr void Load(const float* vB) { this->v = vld1q_f32(vB); } + constexpr void Store(float* vB) const { vst1q_f32(vB, this->v); } + + // fp16 <-> fp32 conversion (fcvtl / fcvtn) is part of the AArch64 + // baseline. The halves move through memcpy because arm_neon.h types + // its fp16 pointers as __fp16, which _Float16 may not alias. + constexpr void Load(const _Float16* vB) { + float16x4_t half; + std::memcpy(&half, vB, sizeof(half)); + this->v = vcvt_f32_f16(half); + } + constexpr void Store(_Float16* vB) const { + float16x4_t half = vcvt_f16_f32(this->v); + std::memcpy(vB, &half, sizeof(half)); + } + + template + constexpr std::array Store() const { + std::array r; + Store(r.data()); + return r; + } + + template + constexpr operator VectorF32() const { + return VectorF32(Permute()>(this->v, this->v)); + } + + constexpr VectorF32 operator+(VectorF32 b) const { return VectorF32(vaddq_f32(this->v, b.v)); } + constexpr VectorF32 operator-(VectorF32 b) const { return VectorF32(vsubq_f32(this->v, b.v)); } + constexpr VectorF32 operator*(VectorF32 b) const { return VectorF32(vmulq_f32(this->v, b.v)); } + constexpr VectorF32 operator/(VectorF32 b) const { return VectorF32(vdivq_f32(this->v, b.v)); } + constexpr void operator+=(VectorF32 b) { this->v = vaddq_f32(this->v, b.v); } + constexpr void operator-=(VectorF32 b) { this->v = vsubq_f32(this->v, b.v); } + constexpr void operator*=(VectorF32 b) { this->v = vmulq_f32(this->v, b.v); } + constexpr void operator/=(VectorF32 b) { this->v = vdivq_f32(this->v, b.v); } + + constexpr VectorF32 operator+(float b) const { return *this + VectorF32(b); } + constexpr VectorF32 operator-(float b) const { return *this - VectorF32(b); } + constexpr VectorF32 operator*(float b) const { return *this * VectorF32(b); } + constexpr VectorF32 operator/(float b) const { return *this / VectorF32(b); } + constexpr void operator+=(float b) { *this += VectorF32(b); } + constexpr void operator-=(float b) { *this -= VectorF32(b); } + constexpr void operator*=(float b) { *this *= VectorF32(b); } + constexpr void operator/=(float b) { *this /= VectorF32(b); } + + constexpr VectorF32 operator-() const { return VectorF32(vnegq_f32(this->v)); } + + // Only the Len*Packing lanes in use take part; padding lanes can hold + // anything. + constexpr bool operator==(VectorF32 b) const { + uint32x4_t equal = vceqq_f32(this->v, b.v); + if constexpr (Len * Packing == 4) { + return vminvq_u32(equal) == ~0u; + } else if constexpr (Len * Packing == 3) { + return vminvq_u32(vcopyq_laneq_u32(equal, 3, equal, 0)) == ~0u; + } else if constexpr (Len * Packing == 2) { + return vgetq_lane_u64(vreinterpretq_u64_u32(equal), 0) == ~0ull; + } else { + return vgetq_lane_u32(equal, 0) == ~0u; + } + } + constexpr bool operator!=(VectorF32 b) const { return !(*this == b); } + + template + constexpr VectorF32 ExtractLo() const { + return VectorF32(Permute()>(this->v, this->v)); + } + + constexpr VectorF32 Cos() const { + return VectorF32(Base::cos_f32x4(this->v)); + } + constexpr VectorF32 Sin() const { + return VectorF32(Base::sin_f32x4(this->v)); + } + constexpr std::tuple, VectorF32> SinCos() const { + float32x4_t s, c; + Base::sincos_f32x4(this->v, s, c); + return { VectorF32(s), VectorF32(c) }; + } + + template values> + constexpr VectorF32 Negate() const { + constexpr std::array mask = Base::template GetNegateMask(); + return VectorF32(vreinterpretq_f32_u32(veorq_u32(vreinterpretq_u32_f32(this->v), vreinterpretq_u32_f32(vld1q_f32(mask.data()))))); + } + + // vfmaq_f32(a, b, c) is a + b*c. + static constexpr VectorF32 MulitplyAdd(VectorF32 a, VectorF32 b, VectorF32 add) { + return VectorF32(vfmaq_f32(add.v, a.v, b.v)); + } + static constexpr VectorF32 MulitplySub(VectorF32 a, VectorF32 b, VectorF32 sub) { + return VectorF32(vfmaq_f32(vnegq_f32(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. The rotations let the padding lane + // repeat y, which is two ext instructions; keeping that lane in place, + // as Shuffle does, takes three. vfmsq_f32(a, b, c) is a - b*c. + constexpr static VectorF32 Cross(VectorF32 a, VectorF32 b) requires(Len == 3) { + constexpr LaneMap yzx{1, 2, 0, 1}; + float32x4_t aYzx = Permute(a.v, a.v); + float32x4_t bYzx = Permute(b.v, b.v); + float32x4_t zxy = vfmsq_f32(vmulq_f32(a.v, bYzx), b.v, aYzx); + return VectorF32(Permute(zxy, zxy)); + } + + template ShuffleValues> + constexpr VectorF32 Shuffle() const { + return VectorF32(Permute()>(this->v, this->v)); + } + + template ShuffleValues> + constexpr static VectorF32 Blend(VectorF32 a, VectorF32 b) { + return VectorF32(Permute()>(a.v, b.v)); + } + + template + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == VectorBase::BatchSize)) + constexpr static auto LengthSq(VectorF32 first, Rest... rest) { + constexpr std::uint8_t N = VectorBase::BatchSize; + std::array, N> args{ first, rest... }; + std::array squares; + for (std::uint8_t i = 0; i < N; ++i) squares[i] = vmulq_f32(args[i].v, args[i].v); + return VectorF32<1, static_cast(Packing * N)>(SumPacks(squares)); + } + + template + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == VectorBase::BatchSize)) + constexpr static auto Length(VectorF32 first, Rest... rest) { + auto sq = LengthSq(first, rest...); + sq.v = vsqrtq_f32(sq.v); + return sq; + } + + // Pairwise dot products across BatchSize pairs: Dot(a0, b0, a1, b1, + // ...) returns one lane per sub-vector pair. + template + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == 2 * VectorBase::BatchSize)) + constexpr static auto Dot(VectorF32 first, Rest... rest) { + constexpr std::uint8_t N = VectorBase::BatchSize; + std::array, 2 * N> args{ first, rest... }; + std::array products; + for (std::uint8_t i = 0; i < N; ++i) products[i] = vmulq_f32(args[2 * i].v, args[2 * i + 1].v); + return VectorF32<1, static_cast(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 + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == VectorBase::BatchSize)) + constexpr static auto Normalize(VectorF32 first, Rest... rest) { + constexpr std::uint8_t N = VectorBase::BatchSize; + std::array, N> args{ first, rest... }; + float32x4_t scale = vdivq_f32(vdupq_n_f32(1.0f), Length(first, rest...).v); + return [&](std::index_sequence) { + return std::array, N>{ + VectorF32(vmulq_f32(args[Is].v, Permute()>(scale, scale)))... + }; + }(std::make_index_sequence{}); + } + + constexpr static VectorF32 Rotate(VectorF32<3, Packing> v, VectorF32<4, Packing> q) requires(Len == 3) { + VectorF32<3, Packing> qv = q; + VectorF32<3, Packing> qw = q.template Shuffle<{{3,3,3,3}}>(); + VectorF32<3, Packing> t = Cross(qv, v) * 2.0f; + return MulitplyAdd(t, qw, v) + Cross(qv, t); + } + + constexpr static VectorF32<3, Packing> RotatePivot(VectorF32<3, Packing> v, VectorF32<4, Packing> q, VectorF32<3, Packing> pivot) requires(Len == 3) { + VectorF32<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 (lanes + // 0..3) and cos (lanes 4..7) registers. + constexpr static VectorF32<4, Packing> QuanternionFromEuler(VectorF32<3, Packing> EulerHalf) requires(Len == 4) { + float32x4_t s, c; + Base::sincos_f32x4(EulerHalf.v, s, c); + float32x4_t left = vmulq_f32( + vmulq_f32( + Permute(s, c), // sx cx cx cx + Permute(s, c)), // cy sy cy cy + Permute(s, c)); // cz cz sz cz + float32x4_t right = vmulq_f32( + Permute(s, c), // cx sx sx sx + Permute(s, c)); // sy cy sy sy + VectorF32<4, Packing> rightZ(Permute(s, c)); // sz sz cz sz + rightZ = rightZ.template Negate<{{true, false, true, false}}>(); + return VectorF32<4, Packing>(vfmaq_f32(left, right, rightZ.v)); + } + }; #else - // Scalar software fallback for non-x86_64 targets. Future arches can swap - // in their own intrinsic implementation by adding an arch-specific branch - // above and gating this one out. + // Scalar software fallback for targets without an intrinsic backend above. + // Future arches can swap in their own intrinsic implementation by adding + // an arch-specific branch above and gating this one out. export template struct VectorF32 : public VectorBase { template @@ -2136,6 +2447,28 @@ namespace Crafter { return sq; } + // Batched dot product over BatchSize (a, b) pairs: Dot(a0, b0, a1, b1, + // ...) returns one lane per pair. MatrixRowMajor relies on it, so + // without it this fallback (every target that has no intrinsic + // backend above) does not compile. Only the first Len lanes of each + // input are summed, as in the other backends. + template + requires((std::is_same_v> && ...) && + (1 + sizeof...(Rest) == 2 * VectorBase::BatchSize)) + constexpr static auto Dot(VectorF32 first, Rest... rest) { + constexpr std::uint8_t N = VectorBase::BatchSize; + VectorF32<1, static_cast(Packing * N)> r; + std::array, 2 * N> args{ first, rest... }; + for (std::uint8_t i = 0; i < N; ++i) + for (std::uint8_t p = 0; p < Packing; ++p) { + float acc = 0.0f; + for (std::uint8_t k = 0; k < Len; ++k) + acc += args[2 * i].v[p * Len + k] * args[2 * i + 1].v[p * Len + k]; + r.v[i * Packing + p] = acc; + } + return r; + } + template requires((std::is_same_v> && ...) && (1 + sizeof...(Rest) == VectorBase::BatchSize)) diff --git a/project.cpp b/project.cpp index e07becf..be528b6 100644 --- a/project.cpp +++ b/project.cpp @@ -46,12 +46,22 @@ extern "C" Configuration CrafterBuildProject(std::span a {"rv64gcv", "generic"}, }}; + // AArch64 has a single NEON register width (128-bit), so the tiers + // differ in VectorF16 only: ARMv8.2-FP16 gives it a native half-precision + // path, the baseline aliases it to VectorF32. + std::array armTiers = {{ + {"armv8.2-a+fp16", "generic"}, + {"armv8-a", "generic"}, + }}; + const std::string_view target = cfg.target; const bool isX86 = target.starts_with("x86_64") || target.starts_with("i686"); const bool isRiscv = target.starts_with("riscv"); + const bool isArm = target.starts_with("aarch64"); for (std::string_view name : { "Vector", "Intersection", "Matrix" }) { if (isX86) cfg.AddMarchVariants(name, mathInterfaces, x86Tiers); if (isRiscv) cfg.AddMarchVariants(name, mathInterfaces, riscvTiers); + if (isArm) cfg.AddMarchVariants(name, mathInterfaces, armTiers); } return cfg;