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 <noreply@anthropic.com>

.gitattributes

@@ -0,0 +1 @@
+* text=auto eol=lf

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<T, Len, Alignment>` with `.x/.y/.z/.w` and `.r/.g/.b/.a` accessors, arithmetic, comparison, and conversion operators. |
 | `:MatrixRowMajor` | `MatrixRowMajor<T, ColumnSize, RowSize, Repeats>` — row-major matrices with packed-batch support via `Repeats`. |
 | `:VectorF32` | `VectorF32<Len, Packing>` — SIMD-specialized 32-bit float vector with `Load`/`Store`, FMA, `Cross`, `Dot`, `Length`, `Normalize`, `Sin`/`Cos`, `Shuffle`, etc. |
 | `:VectorF16` | `VectorF16<Len, Packing>` — 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<T, N>`.
 

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 <std::uint32_t Len, std::uint32_t Packing>
     using VectorF16 = VectorF32<Len, Packing>;
     #endif

interfaces/Crafter.Math-Common.cppm

@@ -8,6 +8,9 @@ module;
 #ifdef __wasm_simd128__
 #include <wasm_simd128.h>
 #endif
+#if defined(__aarch64__) && defined(__ARM_NEON)
+#include <arm_neon.h>
+#endif
 #ifdef __riscv_vector
 #include <riscv_vector.h>
 // 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
+// a Zvfh path lands, and on AArch64 NEON without ARMv8.2-FP16 arithmetic,
-// redefines the same condition since macros do not cross module boundaries.
+// where VectorF32 converts fp16 <-> fp32 on Load/Store. Each translation unit
-#if defined(__FLT16_MAX__) && (!defined(__x86_64) || defined(__AVX512FP16__)) && !defined(__riscv_vector)
+// 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 <std::uint8_t Len, std::uint8_t Packing>
     struct VectorF16;
@@ -52,7 +57,7 @@ namespace Crafter {
 
     template <std::uint8_t Len, std::uint8_t Packing, typename T>
     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 <std::uint8_t L, std::uint8_t P>
         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<std::is_same_v<T, float>, 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
     };
 }

interfaces/Crafter.Math-VectorF16.cppm

@@ -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
+// are on an x86_64 build with AVX512FP16 or an AArch64 NEON build with
-// non-x86_64 target (scalar fallback further below). On x86_64 without
+// ARMv8.2-FP16 (intrinsic impls below), or on a target with neither x86_64 nor
-// AVX512FP16 or anywhere without _Float16, Crafter.Math:Basic aliases
+// NEON (scalar fallback further below). On x86_64 without AVX512FP16, on
-// VectorF16 to VectorF32 instead.
+// 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 {

interfaces/Crafter.Math-VectorF32.cppm

@@ -11,6 +11,9 @@ module;
 #ifdef __riscv_vector
 #include <riscv_vector.h>
 #endif
+#if defined(__aarch64__) && defined(__ARM_NEON)
+#include <arm_neon.h>
+#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 <std::uint8_t Len, std::uint8_t Packing>
+    struct VectorF32 : public VectorBase<Len, Packing, float> {
+        template <std::uint8_t Len2, std::uint8_t Packing2>
+        friend struct VectorF32;
+        using Base = VectorBase<Len, Packing, float>;
+        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<int, 4>;
+
+        template <LaneMap Lanes>
+        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 <std::array<std::uint8_t, Len> 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 <std::array<bool, Len> 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 <std::uint8_t ToLen, std::uint8_t ToPacking>
+        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 <std::uint8_t Index>
+        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<float32x4_t, Base::BatchSize> 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<typename T>
+        constexpr std::array<T, NElems> Store() const {
+            std::array<T, NElems> r;
+            Store(r.data());
+            return r;
+        }
+
+        template <std::uint8_t BLen, std::uint8_t BPacking>
+        constexpr operator VectorF32<BLen, BPacking>() const {
+            return VectorF32<BLen, BPacking>(Permute<GetResizeLanes<BLen, BPacking>()>(this->v, this->v));
+        }
+
+        constexpr VectorF32<Len, Packing> operator+(VectorF32<Len, Packing> b) const { return VectorF32<Len, Packing>(vaddq_f32(this->v, b.v)); }
+        constexpr VectorF32<Len, Packing> operator-(VectorF32<Len, Packing> b) const { return VectorF32<Len, Packing>(vsubq_f32(this->v, b.v)); }
+        constexpr VectorF32<Len, Packing> operator*(VectorF32<Len, Packing> b) const { return VectorF32<Len, Packing>(vmulq_f32(this->v, b.v)); }
+        constexpr VectorF32<Len, Packing> operator/(VectorF32<Len, Packing> b) const { return VectorF32<Len, Packing>(vdivq_f32(this->v, b.v)); }
+        constexpr void operator+=(VectorF32<Len, Packing> b) { this->v = vaddq_f32(this->v, b.v); }
+        constexpr void operator-=(VectorF32<Len, Packing> b) { this->v = vsubq_f32(this->v, b.v); }
+        constexpr void operator*=(VectorF32<Len, Packing> b) { this->v = vmulq_f32(this->v, b.v); }
+        constexpr void operator/=(VectorF32<Len, Packing> b) { this->v = vdivq_f32(this->v, b.v); }
+
+        constexpr VectorF32<Len, Packing> operator+(float b) const { return *this + VectorF32<Len, Packing>(b); }
+        constexpr VectorF32<Len, Packing> operator-(float b) const { return *this - VectorF32<Len, Packing>(b); }
+        constexpr VectorF32<Len, Packing> operator*(float b) const { return *this * VectorF32<Len, Packing>(b); }
+        constexpr VectorF32<Len, Packing> operator/(float b) const { return *this / VectorF32<Len, Packing>(b); }
+        constexpr void operator+=(float b) { *this += VectorF32<Len, Packing>(b); }
+        constexpr void operator-=(float b) { *this -= VectorF32<Len, Packing>(b); }
+        constexpr void operator*=(float b) { *this *= VectorF32<Len, Packing>(b); }
+        constexpr void operator/=(float b) { *this /= VectorF32<Len, Packing>(b); }
+
+        constexpr VectorF32<Len, Packing> operator-() const { return VectorF32<Len, Packing>(vnegq_f32(this->v)); }
+
+        // Only the Len*Packing lanes in use take part; padding lanes can hold
+        // anything.
+        constexpr bool operator==(VectorF32<Len, Packing> 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<Len, Packing> b) const { return !(*this == b); }
+
+        template<std::uint32_t ExtractLen>
+        constexpr VectorF32<ExtractLen, Packing> ExtractLo() const {
+            return VectorF32<ExtractLen, Packing>(Permute<GetResizeLanes<ExtractLen, Packing>()>(this->v, this->v));
+        }
+
+        constexpr VectorF32<Len, Packing> Cos() const {
+            return VectorF32<Len, Packing>(Base::cos_f32x4(this->v));
+        }
+        constexpr VectorF32<Len, Packing> Sin() const {
+            return VectorF32<Len, Packing>(Base::sin_f32x4(this->v));
+        }
+        constexpr std::tuple<VectorF32<Len, Packing>, VectorF32<Len, Packing>> SinCos() const {
+            float32x4_t s, c;
+            Base::sincos_f32x4(this->v, s, c);
+            return { VectorF32<Len, Packing>(s), VectorF32<Len, Packing>(c) };
+        }
+
+        template <std::array<bool, Len> values>
+        constexpr VectorF32<Len, Packing> Negate() const {
+            constexpr std::array<float, NElems> mask = Base::template GetNegateMask<values>();
+            return VectorF32<Len, Packing>(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<Len, Packing> MulitplyAdd(VectorF32<Len, Packing> a, VectorF32<Len, Packing> b, VectorF32<Len, Packing> add) {
+            return VectorF32<Len, Packing>(vfmaq_f32(add.v, a.v, b.v));
+        }
+        static constexpr VectorF32<Len, Packing> MulitplySub(VectorF32<Len, Packing> a, VectorF32<Len, Packing> b, VectorF32<Len, Packing> sub) {
+            return VectorF32<Len, Packing>(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<Len, Packing> Cross(VectorF32<Len, Packing> a, VectorF32<Len, Packing> b) requires(Len == 3) {
+            constexpr LaneMap yzx{1, 2, 0, 1};
+            float32x4_t aYzx = Permute<yzx>(a.v, a.v);
+            float32x4_t bYzx = Permute<yzx>(b.v, b.v);
+            float32x4_t zxy = vfmsq_f32(vmulq_f32(a.v, bYzx), b.v, aYzx);
+            return VectorF32<Len, Packing>(Permute<yzx>(zxy, zxy));
+        }
+
+        template <const std::array<std::uint8_t, Len> ShuffleValues>
+        constexpr VectorF32<Len, Packing> Shuffle() const {
+            return VectorF32<Len, Packing>(Permute<GetShuffleLanes<ShuffleValues>()>(this->v, this->v));
+        }
+
+        template <std::array<bool, Len> ShuffleValues>
+        constexpr static VectorF32<Len, Packing> Blend(VectorF32<Len, Packing> a, VectorF32<Len, Packing> b) {
+            return VectorF32<Len, Packing>(Permute<GetBlendLanes<ShuffleValues>()>(a.v, b.v));
+        }
+
+        template<typename... Rest>
+            requires((std::is_same_v<Rest, VectorF32<Len, Packing>> && ...) &&
+                     (1 + sizeof...(Rest) == VectorBase<Len, Packing, float>::BatchSize))
+        constexpr static auto LengthSq(VectorF32<Len, Packing> first, Rest... rest) {
+            constexpr std::uint8_t N = VectorBase<Len, Packing, float>::BatchSize;
+            std::array<VectorF32<Len, Packing>, N> args{ first, rest... };
+            std::array<float32x4_t, N> 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<std::uint8_t>(Packing * N)>(SumPacks(squares));
+        }
+
+        template<typename... Rest>
+            requires((std::is_same_v<Rest, VectorF32<Len, Packing>> && ...) &&
+                     (1 + sizeof...(Rest) == VectorBase<Len, Packing, float>::BatchSize))
+        constexpr static auto Length(VectorF32<Len, Packing> 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<typename... Rest>
+            requires((std::is_same_v<Rest, VectorF32<Len, Packing>> && ...) &&
+                     (1 + sizeof...(Rest) == 2 * VectorBase<Len, Packing, float>::BatchSize))
+        constexpr static auto Dot(VectorF32<Len, Packing> first, Rest... rest) {
+            constexpr std::uint8_t N = VectorBase<Len, Packing, float>::BatchSize;
+            std::array<VectorF32<Len, Packing>, 2 * N> args{ first, rest... };
+            std::array<float32x4_t, N> 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<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, VectorF32<Len, Packing>> && ...) &&
+                     (1 + sizeof...(Rest) == VectorBase<Len, Packing, float>::BatchSize))
+        constexpr static auto Normalize(VectorF32<Len, Packing> first, Rest... rest) {
+            constexpr std::uint8_t N = VectorBase<Len, Packing, float>::BatchSize;
+            std::array<VectorF32<Len, Packing>, N> args{ first, rest... };
+            float32x4_t scale = vdivq_f32(vdupq_n_f32(1.0f), Length(first, rest...).v);
+            return [&]<std::size_t... Is>(std::index_sequence<Is...>) {
+                return std::array<VectorF32<Len, Packing>, N>{
+                    VectorF32<Len, Packing>(vmulq_f32(args[Is].v, Permute<GetSpreadLanes<Is>()>(scale, scale)))...
+                };
+            }(std::make_index_sequence<N>{});
+        }
+
+        constexpr static VectorF32<Len, Packing> 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<LaneMap{0, 4, 4, 4}>(s, c),     // sx cx cx cx
+                    Permute<LaneMap{5, 1, 5, 5}>(s, c)),    // cy sy cy cy
+                Permute<LaneMap{6, 6, 2, 6}>(s, c));        // cz cz sz cz
+            float32x4_t right = vmulq_f32(
+                Permute<LaneMap{4, 0, 0, 0}>(s, c),         // cx sx sx sx
+                Permute<LaneMap{1, 5, 1, 1}>(s, c));        // sy cy sy sy
+            VectorF32<4, Packing> rightZ(Permute<LaneMap{2, 2, 6, 2}>(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
+    // Scalar software fallback for targets without an intrinsic backend above.
-    // in their own intrinsic implementation by adding an arch-specific branch
+    // Future arches can swap in their own intrinsic implementation by adding
-    // above and gating this one out.
+    // an arch-specific branch above and gating this one out.
     export template <std::uint8_t Len, std::uint8_t Packing>
     struct VectorF32 : public VectorBase<Len, Packing, float> {
         template <std::uint8_t Len2, std::uint8_t Packing2>
@@ -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<typename... Rest>
+            requires((std::is_same_v<Rest, VectorF32<Len, Packing>> && ...) &&
+                     (1 + sizeof...(Rest) == 2 * VectorBase<Len, Packing, float>::BatchSize))
+        constexpr static auto Dot(VectorF32<Len, Packing> first, Rest... rest) {
+            constexpr std::uint8_t N = VectorBase<Len, Packing, float>::BatchSize;
+            VectorF32<1, static_cast<std::uint8_t>(Packing * N)> r;
+            std::array<VectorF32<Len, Packing>, 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<typename... Rest>
             requires((std::is_same_v<Rest, VectorF32<Len, Packing>> && ...) &&
                      (1 + sizeof...(Rest) == VectorBase<Len, Packing, float>::BatchSize))

project.cpp

@@ -46,12 +46,22 @@ extern "C" Configuration CrafterBuildProject(std::span<const std::string_view> 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<MarchTier, 2> 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;