From cfc54bf1f74dcab84db7b491cc5b703409b7f1d1 Mon Sep 17 00:00:00 2001 From: Miha Kralj Date: Wed, 10 Dec 2025 21:12:14 -0500 Subject: [PATCH] Enhance SIMD support in moving average implementations; add AVX512 and ARM64 optimizations for HMA, SMA, and WMA --- lib/core/simd/SimdExtensions.cs | 73 ++------- lib/trends/hma/Hma.cs | 39 +++-- lib/trends/sma/Sma.cs | 165 ++++++++++++++++++- lib/trends/wma/Wma.cs | 276 ++++++++++++++++++++++++++++++++ lib/trends/wma/Wma.md | 1 + 5 files changed, 476 insertions(+), 78 deletions(-) diff --git a/lib/core/simd/SimdExtensions.cs b/lib/core/simd/SimdExtensions.cs index 088a214c..fe74a931 100644 --- a/lib/core/simd/SimdExtensions.cs +++ b/lib/core/simd/SimdExtensions.cs @@ -383,7 +383,7 @@ public static class SimdExtensions /// /// Calculates the dot product of two spans using SIMD intrinsics. - /// Supports AVX512, AVX2, SSE2, and NEON. + /// Supports AVX512, AVX2, and NEON (ARM64). /// [MethodImpl(MethodImplOptions.AggressiveInlining)] public static double DotProduct(this ReadOnlySpan a, ReadOnlySpan b) @@ -410,26 +410,25 @@ public static class SimdExtensions if (Avx2.IsSupported) return DotProductAvx2(a, b); - if (Sse2.IsSupported) - return DotProductSse2(a, b); - if (AdvSimd.Arm64.IsSupported) return DotProductNeon(a, b); - double s = 0; + double s1 = 0, s2 = 0, s3 = 0, s4 = 0; ref double ar = ref MemoryMarshal.GetReference(a); ref double br = ref MemoryMarshal.GetReference(b); int i = 0; - // Unroll scalar loop + // Unroll scalar loop with 4 accumulators to break dependency chains for (; i <= len - 4; i += 4) { - s += Unsafe.Add(ref ar, i) * Unsafe.Add(ref br, i); - s += Unsafe.Add(ref ar, i + 1) * Unsafe.Add(ref br, i + 1); - s += Unsafe.Add(ref ar, i + 2) * Unsafe.Add(ref br, i + 2); - s += Unsafe.Add(ref ar, i + 3) * Unsafe.Add(ref br, i + 3); + s1 += Unsafe.Add(ref ar, i) * Unsafe.Add(ref br, i); + s2 += Unsafe.Add(ref ar, i + 1) * Unsafe.Add(ref br, i + 1); + s3 += Unsafe.Add(ref ar, i + 2) * Unsafe.Add(ref br, i + 2); + s4 += Unsafe.Add(ref ar, i + 3) * Unsafe.Add(ref br, i + 3); } + double s = s1 + s2 + s3 + s4; + for (; i < len; i++) { s += Unsafe.Add(ref ar, i) * Unsafe.Add(ref br, i); @@ -647,58 +646,4 @@ public static class SimdExtensions return sum; } - [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static double DotProductSse2(ReadOnlySpan a, ReadOnlySpan b) - { - int len = a.Length; - int i = 0; - Vector128 vSum = Vector128.Zero; - Vector128 vSum2 = Vector128.Zero; - - ref double aRef = ref MemoryMarshal.GetReference(a); - ref double bRef = ref MemoryMarshal.GetReference(b); - - // Process 4 doubles at a time using 2 accumulators - for (; i <= len - 4; i += 4) - { - var va1 = Vector128.LoadUnsafe(ref Unsafe.Add(ref aRef, i)); - var vb1 = Vector128.LoadUnsafe(ref Unsafe.Add(ref bRef, i)); - var va2 = Vector128.LoadUnsafe(ref Unsafe.Add(ref aRef, i + 2)); - var vb2 = Vector128.LoadUnsafe(ref Unsafe.Add(ref bRef, i + 2)); - - if (Fma.IsSupported) - { - vSum = Fma.MultiplyAdd(va1, vb1, vSum); - vSum2 = Fma.MultiplyAdd(va2, vb2, vSum2); - } - else - { - vSum = Sse2.Add(vSum, Sse2.Multiply(va1, vb1)); - vSum2 = Sse2.Add(vSum2, Sse2.Multiply(va2, vb2)); - } - } - - // Process remaining 2 doubles if available - if (i <= len - 2) - { - var va = Vector128.LoadUnsafe(ref Unsafe.Add(ref aRef, i)); - var vb = Vector128.LoadUnsafe(ref Unsafe.Add(ref bRef, i)); - - vSum = Fma.IsSupported - ? Fma.MultiplyAdd(va, vb, vSum) - : Sse2.Add(vSum, Sse2.Multiply(va, vb)); - i += 2; - } - - vSum = Sse2.Add(vSum, vSum2); - double sum = vSum.GetElement(0) + vSum.GetElement(1); - - // Scalar remainder (0-1 elements) - for (; i < len; i++) - { - sum += Unsafe.Add(ref aRef, i) * Unsafe.Add(ref bRef, i); - } - - return sum; - } } diff --git a/lib/trends/hma/Hma.cs b/lib/trends/hma/Hma.cs index ebb4bcc4..8edf38f8 100644 --- a/lib/trends/hma/Hma.cs +++ b/lib/trends/hma/Hma.cs @@ -2,6 +2,8 @@ using System; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Runtime.Intrinsics; +using System.Runtime.Intrinsics.Arm; +using System.Runtime.Intrinsics.X86; namespace QuanTAlib; @@ -163,24 +165,43 @@ public sealed class Hma : ITValuePublisher int len = halfWma.Length; int i = 0; - if (Vector256.IsHardwareAccelerated && len >= Vector256.Count) + ref double halfRef = ref MemoryMarshal.GetReference(halfWma); + ref double fullRef = ref MemoryMarshal.GetReference(fullWma); + ref double outRef = ref MemoryMarshal.GetReference(output); + + if (Avx512F.IsSupported && len >= Vector512.Count) + { + var vTwo = Vector512.Create(2.0); + for (; i <= len - Vector512.Count; i += Vector512.Count) + { + var vHalf = Vector512.LoadUnsafe(ref Unsafe.Add(ref halfRef, i)); + var vFull = Vector512.LoadUnsafe(ref Unsafe.Add(ref fullRef, i)); + var vResult = Avx512F.Subtract(Avx512F.Multiply(vHalf, vTwo), vFull); + Vector512.StoreUnsafe(vResult, ref Unsafe.Add(ref outRef, i)); + } + } + else if (Avx2.IsSupported && len >= Vector256.Count) { var vTwo = Vector256.Create(2.0); - ref double halfRef = ref MemoryMarshal.GetReference(halfWma); - ref double fullRef = ref MemoryMarshal.GetReference(fullWma); - ref double outRef = ref MemoryMarshal.GetReference(output); - for (; i <= len - Vector256.Count; i += Vector256.Count) { var vHalf = Vector256.LoadUnsafe(ref Unsafe.Add(ref halfRef, i)); var vFull = Vector256.LoadUnsafe(ref Unsafe.Add(ref fullRef, i)); - - // vResult = 2 * half - full - var vResult = (vHalf * vTwo) - vFull; - + var vResult = Avx.Subtract(Avx.Multiply(vHalf, vTwo), vFull); Vector256.StoreUnsafe(vResult, ref Unsafe.Add(ref outRef, i)); } } + else if (AdvSimd.Arm64.IsSupported && len >= Vector128.Count) + { + var vTwo = Vector128.Create(2.0); + for (; i <= len - Vector128.Count; i += Vector128.Count) + { + var vHalf = Vector128.LoadUnsafe(ref Unsafe.Add(ref halfRef, i)); + var vFull = Vector128.LoadUnsafe(ref Unsafe.Add(ref fullRef, i)); + var vResult = AdvSimd.Arm64.Subtract(AdvSimd.Arm64.Multiply(vHalf, vTwo), vFull); + Vector128.StoreUnsafe(vResult, ref Unsafe.Add(ref outRef, i)); + } + } for (; i < len; i++) { diff --git a/lib/trends/sma/Sma.cs b/lib/trends/sma/Sma.cs index e442de07..86113934 100644 --- a/lib/trends/sma/Sma.cs +++ b/lib/trends/sma/Sma.cs @@ -3,6 +3,7 @@ using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Runtime.Intrinsics; +using System.Runtime.Intrinsics.Arm; using System.Runtime.Intrinsics.X86; namespace QuanTAlib; @@ -223,12 +224,27 @@ public sealed class Sma : ITValuePublisher if (len == 0) return; // Try SIMD path for large, clean datasets - // Requirements: AVX2 support, large enough dataset, no NaN values + // Requirements: SIMD support, large enough dataset, no NaN values const int SimdThreshold = 256; - if (Avx2.IsSupported && len >= SimdThreshold && !source.ContainsNonFinite()) + if (len >= SimdThreshold && !source.ContainsNonFinite()) { - CalculateSimdCore(source, output, period); - return; + if (Avx512F.IsSupported) + { + CalculateAvx512Core(source, output, period); + return; + } + + if (Avx2.IsSupported) + { + CalculateAvx2Core(source, output, period); + return; + } + + if (AdvSimd.Arm64.IsSupported) + { + CalculateNeonCore(source, output, period); + return; + } } // Scalar path with NaN handling @@ -297,7 +313,79 @@ public sealed class Sma : ITValuePublisher } [MethodImpl(MethodImplOptions.AggressiveOptimization)] - private static void CalculateSimdCore(ReadOnlySpan source, Span output, int period) + private static void CalculateAvx512Core(ReadOnlySpan source, Span output, int period) + { + int len = source.Length; + const int VectorWidth = 8; + + ref double srcRef = ref MemoryMarshal.GetReference(source); + ref double outRef = ref MemoryMarshal.GetReference(output); + + double invPeriod = 1.0 / period; + + int warmupEnd = Math.Min(period, len); + double sum = 0; + for (int i = 0; i < warmupEnd; i++) + { + sum += Unsafe.Add(ref srcRef, i); + Unsafe.Add(ref outRef, i) = sum / (i + 1); + } + + if (len <= period) + return; + + var vInvPeriod = Vector512.Create(invPeriod); + int simdEnd = period + ((len - period) / VectorWidth) * VectorWidth; + int tickCount = 0; + + for (int i = period; i < simdEnd; i += VectorWidth) + { + var vNew = Vector512.LoadUnsafe(ref Unsafe.Add(ref srcRef, i)); + var vOld = Vector512.LoadUnsafe(ref Unsafe.Add(ref srcRef, i - period)); + + var vDelta = Avx512F.Subtract(vNew, vOld); + + // Prefix sum of Delta + var vShift1 = Vector512.Create(0.0, vDelta.GetElement(0), vDelta.GetElement(1), vDelta.GetElement(2), vDelta.GetElement(3), vDelta.GetElement(4), vDelta.GetElement(5), vDelta.GetElement(6)); + var vP1 = Avx512F.Add(vDelta, vShift1); + + var vShift2 = Vector512.Create(0.0, 0.0, vP1.GetElement(0), vP1.GetElement(1), vP1.GetElement(2), vP1.GetElement(3), vP1.GetElement(4), vP1.GetElement(5)); + var vP2 = Avx512F.Add(vP1, vShift2); + + var vShift4 = Vector512.Create(0.0, 0.0, 0.0, 0.0, vP2.GetElement(0), vP2.GetElement(1), vP2.GetElement(2), vP2.GetElement(3)); + var vP4 = Avx512F.Add(vP2, vShift4); + + var vSumPrev = Vector512.Create(sum); + var vSums = Avx512F.Add(vSumPrev, vP4); + + var vResult = Avx512F.Multiply(vSums, vInvPeriod); + Vector512.StoreUnsafe(vResult, ref Unsafe.Add(ref outRef, i)); + + sum = vSums.GetElement(7); + + tickCount += VectorWidth; + if (tickCount >= ResyncInterval) + { + tickCount = 0; + int lastIdx = i + VectorWidth - 1; + double recalcSum = 0; + for (int k = 0; k < period; k++) + { + recalcSum += Unsafe.Add(ref srcRef, lastIdx - k); + } + sum = recalcSum; + } + } + + for (int i = simdEnd; i < len; i++) + { + sum = sum - Unsafe.Add(ref srcRef, i - period) + Unsafe.Add(ref srcRef, i); + Unsafe.Add(ref outRef, i) = sum * invPeriod; + } + } + + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void CalculateAvx2Core(ReadOnlySpan source, Span output, int period) { int len = source.Length; const int VectorWidth = 4; @@ -367,6 +455,73 @@ public sealed class Sma : ITValuePublisher } } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void CalculateNeonCore(ReadOnlySpan source, Span output, int period) + { + int len = source.Length; + const int VectorWidth = 2; + + ref double srcRef = ref MemoryMarshal.GetReference(source); + ref double outRef = ref MemoryMarshal.GetReference(output); + + double invPeriod = 1.0 / period; + + int warmupEnd = Math.Min(period, len); + double sum = 0; + for (int i = 0; i < warmupEnd; i++) + { + sum += Unsafe.Add(ref srcRef, i); + Unsafe.Add(ref outRef, i) = sum / (i + 1); + } + + if (len <= period) + return; + + var vInvPeriod = Vector128.Create(invPeriod); + int simdEnd = period + ((len - period) / VectorWidth) * VectorWidth; + int tickCount = 0; + + for (int i = period; i < simdEnd; i += VectorWidth) + { + var vNew = Vector128.LoadUnsafe(ref Unsafe.Add(ref srcRef, i)); + var vOld = Vector128.LoadUnsafe(ref Unsafe.Add(ref srcRef, i - period)); + + var vDelta = AdvSimd.Arm64.Subtract(vNew, vOld); + + // Prefix sum of Delta: [d0, d0+d1] + double d0 = vDelta.GetElement(0); + double d1 = vDelta.GetElement(1); + double ps0 = sum + d0; + double ps1 = ps0 + d1; + + var vSums = Vector128.Create(ps0, ps1); + + var vResult = AdvSimd.Arm64.Multiply(vSums, vInvPeriod); + Vector128.StoreUnsafe(vResult, ref Unsafe.Add(ref outRef, i)); + + sum = ps1; + + tickCount += VectorWidth; + if (tickCount >= ResyncInterval) + { + tickCount = 0; + int lastIdx = i + VectorWidth - 1; + double recalcSum = 0; + for (int k = 0; k < period; k++) + { + recalcSum += Unsafe.Add(ref srcRef, lastIdx - k); + } + sum = recalcSum; + } + } + + for (int i = simdEnd; i < len; i++) + { + sum = sum - Unsafe.Add(ref srcRef, i - period) + Unsafe.Add(ref srcRef, i); + Unsafe.Add(ref outRef, i) = sum * invPeriod; + } + } + /// /// Resets the SMA state. /// diff --git a/lib/trends/wma/Wma.cs b/lib/trends/wma/Wma.cs index c5b14e8d..395567ff 100644 --- a/lib/trends/wma/Wma.cs +++ b/lib/trends/wma/Wma.cs @@ -3,6 +3,7 @@ using System.Numerics; using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using System.Runtime.Intrinsics; +using System.Runtime.Intrinsics.Arm; using System.Runtime.Intrinsics.X86; namespace QuanTAlib; @@ -213,12 +214,24 @@ public sealed class Wma : ITValuePublisher if (len == 0) return; const int SimdThreshold = 256; + if (Avx512F.IsSupported && len >= SimdThreshold && !source.ContainsNonFinite()) + { + CalculateAvx512Core(source, output, period); + return; + } + if (Avx2.IsSupported && len >= SimdThreshold && !source.ContainsNonFinite()) { CalculateSimdCore(source, output, period); return; } + if (AdvSimd.Arm64.IsSupported && len >= SimdThreshold && !source.ContainsNonFinite()) + { + CalculateNeonCore(source, output, period); + return; + } + CalculateScalarCore(source, output, period); } @@ -295,6 +308,121 @@ public sealed class Wma : ITValuePublisher } } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void CalculateAvx512Core(ReadOnlySpan source, Span output, int period) + { + int len = source.Length; + const int VectorWidth = 8; + + ref double srcRef = ref MemoryMarshal.GetReference(source); + ref double outRef = ref MemoryMarshal.GetReference(output); + + double divisor = (double)period * (period + 1) * 0.5; + double invDivisor = 1.0 / divisor; + + int warmupEnd = Math.Min(period, len); + double sum = 0; + double wsum = 0; + for (int i = 0; i < warmupEnd; i++) + { + double val = Unsafe.Add(ref srcRef, i); + sum += val; + wsum += (i + 1) * val; + double currentDivisor = (double)(i + 1) * (i + 2) * 0.5; + Unsafe.Add(ref outRef, i) = wsum / currentDivisor; + } + + if (len <= period) + return; + + var vInvDivisor = Vector512.Create(invDivisor); + var vPeriod = Vector512.Create((double)period); + int simdEnd = period + ((len - period) / VectorWidth) * VectorWidth; + + var vSumState = Vector512.Create(sum); + var vWsumState = Vector512.Create(wsum); + + int idx = period; + while (idx < simdEnd) + { + int nextSync = Math.Min(simdEnd, idx + ResyncInterval); + + for (; idx < nextSync; idx += VectorWidth) + { + var vNew = Vector512.LoadUnsafe(ref Unsafe.Add(ref srcRef, idx)); + var vOld = Vector512.LoadUnsafe(ref Unsafe.Add(ref srcRef, idx - period)); + + var vDeltaS = Avx512F.Subtract(vNew, vOld); + + // Prefix sum of DeltaS + var vShiftS1 = Vector512.Create(0.0, vDeltaS.GetElement(0), vDeltaS.GetElement(1), vDeltaS.GetElement(2), vDeltaS.GetElement(3), vDeltaS.GetElement(4), vDeltaS.GetElement(5), vDeltaS.GetElement(6)); + var vPS1 = Avx512F.Add(vDeltaS, vShiftS1); + + var vShiftS2 = Vector512.Create(0.0, 0.0, vPS1.GetElement(0), vPS1.GetElement(1), vPS1.GetElement(2), vPS1.GetElement(3), vPS1.GetElement(4), vPS1.GetElement(5)); + var vPS2 = Avx512F.Add(vPS1, vShiftS2); + + var vShiftS4 = Vector512.Create(0.0, 0.0, 0.0, 0.0, vPS2.GetElement(0), vPS2.GetElement(1), vPS2.GetElement(2), vPS2.GetElement(3)); + var vPS4 = Avx512F.Add(vPS2, vShiftS4); + + var vSums = Avx512F.Add(vSumState, vPS4); + + // Calculate Wsum update + var vSumsShifted = Avx512F.Subtract(vSums, vDeltaS); + var vU = Avx512F.FusedMultiplySubtract(vPeriod, vNew, vSumsShifted); + + // Prefix sum of vU + var vShiftW1 = Vector512.Create(0.0, vU.GetElement(0), vU.GetElement(1), vU.GetElement(2), vU.GetElement(3), vU.GetElement(4), vU.GetElement(5), vU.GetElement(6)); + var vPW1 = Avx512F.Add(vU, vShiftW1); + + var vShiftW2 = Vector512.Create(0.0, 0.0, vPW1.GetElement(0), vPW1.GetElement(1), vPW1.GetElement(2), vPW1.GetElement(3), vPW1.GetElement(4), vPW1.GetElement(5)); + var vPW2 = Avx512F.Add(vPW1, vShiftW2); + + var vShiftW4 = Vector512.Create(0.0, 0.0, 0.0, 0.0, vPW2.GetElement(0), vPW2.GetElement(1), vPW2.GetElement(2), vPW2.GetElement(3)); + var vPW4 = Avx512F.Add(vPW2, vShiftW4); + + var vWsums = Avx512F.Add(vWsumState, vPW4); + + var vResult = Avx512F.Multiply(vWsums, vInvDivisor); + Vector512.StoreUnsafe(vResult, ref Unsafe.Add(ref outRef, idx)); + + // Update state for next iteration + vSumState = Vector512.Create(vSums.GetElement(7)); + vWsumState = Vector512.Create(vWsums.GetElement(7)); + } + + if (idx < len) + { + int lastIdx = idx - 1; + double recalcSum = 0; + double recalcWsum = 0; + for (int k = 0; k < period; k++) + { + double val = Unsafe.Add(ref srcRef, lastIdx - k); + recalcSum += val; + recalcWsum += (period - k) * val; + } + sum = recalcSum; + wsum = recalcWsum; + + vSumState = Vector512.Create(sum); + vWsumState = Vector512.Create(wsum); + } + } + + sum = vSumState.GetElement(0); + wsum = vWsumState.GetElement(0); + + for (; idx < len; idx++) + { + double val = Unsafe.Add(ref srcRef, idx); + double oldSum = sum; + double oldest = Unsafe.Add(ref srcRef, idx - period); + sum = sum - oldest + val; + wsum = wsum - oldSum + (period * val); + Unsafe.Add(ref outRef, idx) = wsum * invDivisor; + } + } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] private static void CalculateSimdCore(ReadOnlySpan source, Span output, int period) { @@ -475,6 +603,154 @@ public sealed class Wma : ITValuePublisher } } + [MethodImpl(MethodImplOptions.AggressiveOptimization)] + private static void CalculateNeonCore(ReadOnlySpan source, Span output, int period) + { + int len = source.Length; + const int VectorWidth = 2; + + ref double srcRef = ref MemoryMarshal.GetReference(source); + ref double outRef = ref MemoryMarshal.GetReference(output); + + double divisor = (double)period * (period + 1) * 0.5; + double invDivisor = 1.0 / divisor; + + int warmupEnd = Math.Min(period, len); + double sum = 0; + double wsum = 0; + for (int i = 0; i < warmupEnd; i++) + { + double val = Unsafe.Add(ref srcRef, i); + sum += val; + wsum += (i + 1) * val; + double currentDivisor = (double)(i + 1) * (i + 2) * 0.5; + Unsafe.Add(ref outRef, i) = wsum / currentDivisor; + } + + if (len <= period) + return; + + var vInvDivisor = Vector128.Create(invDivisor); + int simdEnd = period + ((len - period) / VectorWidth) * VectorWidth; + + double sumState = sum; + double wsumState = wsum; + + int idx = period; + while (idx < simdEnd) + { + int nextSync = Math.Min(simdEnd, idx + ResyncInterval); + + // Unrolled loop: process 4 elements (2 vectors) at a time + int unrolledSync = nextSync - (2 * VectorWidth); + for (; idx <= unrolledSync; idx += 2 * VectorWidth) + { + var vNew1 = Vector128.LoadUnsafe(ref Unsafe.Add(ref srcRef, idx)); + var vOld1 = Vector128.LoadUnsafe(ref Unsafe.Add(ref srcRef, idx - period)); + var vNew2 = Vector128.LoadUnsafe(ref Unsafe.Add(ref srcRef, idx + VectorWidth)); + var vOld2 = Vector128.LoadUnsafe(ref Unsafe.Add(ref srcRef, idx + VectorWidth - period)); + + var vDeltaS1 = AdvSimd.Arm64.Subtract(vNew1, vOld1); + var vDeltaS2 = AdvSimd.Arm64.Subtract(vNew2, vOld2); + + // Prefix sum for first vector: [d0, d0+d1] + double d1_0 = vDeltaS1.GetElement(0); + double d1_1 = vDeltaS1.GetElement(1); + double ps1_0 = sumState + d1_0; + double ps1_1 = ps1_0 + d1_1; + + // Prefix sum for second vector + double d2_0 = vDeltaS2.GetElement(0); + double d2_1 = vDeltaS2.GetElement(1); + double ps2_0 = ps1_1 + d2_0; + double ps2_1 = ps2_0 + d2_1; + + // Calculate Wsum update: W_new = W_old - S_prev + n*new + // For element i: u_i = period * new_i - S_(i-1) + double u1_0 = period * vNew1.GetElement(0) - sumState; + double u1_1 = period * vNew1.GetElement(1) - ps1_0; + double u2_0 = period * vNew2.GetElement(0) - ps1_1; + double u2_1 = period * vNew2.GetElement(1) - ps2_0; + + // Prefix sum of U values + double pw1_0 = wsumState + u1_0; + double pw1_1 = pw1_0 + u1_1; + double pw2_0 = pw1_1 + u2_0; + double pw2_1 = pw2_0 + u2_1; + + var vWsums1 = Vector128.Create(pw1_0, pw1_1); + var vWsums2 = Vector128.Create(pw2_0, pw2_1); + + var vResult1 = AdvSimd.Arm64.Multiply(vWsums1, vInvDivisor); + var vResult2 = AdvSimd.Arm64.Multiply(vWsums2, vInvDivisor); + + Vector128.StoreUnsafe(vResult1, ref Unsafe.Add(ref outRef, idx)); + Vector128.StoreUnsafe(vResult2, ref Unsafe.Add(ref outRef, idx + VectorWidth)); + + sumState = ps2_1; + wsumState = pw2_1; + } + + // Process remaining pairs + for (; idx < nextSync; idx += VectorWidth) + { + var vNew = Vector128.LoadUnsafe(ref Unsafe.Add(ref srcRef, idx)); + var vOld = Vector128.LoadUnsafe(ref Unsafe.Add(ref srcRef, idx - period)); + + var vDeltaS = AdvSimd.Arm64.Subtract(vNew, vOld); + + double d0 = vDeltaS.GetElement(0); + double d1 = vDeltaS.GetElement(1); + double ps0 = sumState + d0; + double ps1 = ps0 + d1; + + double u0 = period * vNew.GetElement(0) - sumState; + double u1 = period * vNew.GetElement(1) - ps0; + + double pw0 = wsumState + u0; + double pw1 = pw0 + u1; + + var vWsums = Vector128.Create(pw0, pw1); + var vResult = AdvSimd.Arm64.Multiply(vWsums, vInvDivisor); + + Vector128.StoreUnsafe(vResult, ref Unsafe.Add(ref outRef, idx)); + + sumState = ps1; + wsumState = pw1; + } + + // Resync to prevent floating-point drift + if (idx < len) + { + int lastIdx = idx - 1; + double recalcSum = 0; + double recalcWsum = 0; + for (int k = 0; k < period; k++) + { + double val = Unsafe.Add(ref srcRef, lastIdx - k); + recalcSum += val; + recalcWsum += (period - k) * val; + } + sumState = recalcSum; + wsumState = recalcWsum; + } + } + + sum = sumState; + wsum = wsumState; + + // Scalar tail + for (; idx < len; idx++) + { + double val = Unsafe.Add(ref srcRef, idx); + double oldSum = sum; + double oldest = Unsafe.Add(ref srcRef, idx - period); + sum = sum - oldest + val; + wsum = wsum - oldSum + (period * val); + Unsafe.Add(ref outRef, idx) = wsum * invDivisor; + } + } + public void Reset() { _buffer.Clear(); diff --git a/lib/trends/wma/Wma.md b/lib/trends/wma/Wma.md index e180cca0..6e5e8a4c 100644 --- a/lib/trends/wma/Wma.md +++ b/lib/trends/wma/Wma.md @@ -211,6 +211,7 @@ The implementation uses: * **Scalar state save/restore** for O(1) bar correction * **Pinned memory** in RingBuffer for cache-friendly access * **CollectionsMarshal.SetCount** for zero-allocation batch processing +* **SIMD Acceleration** (AVX512/AVX2) for high-performance batch processing ## Interpretation Details