From f26c42e9262be4958f5149c74324908144e69829 Mon Sep 17 00:00:00 2001 From: Nobuaki Tanaka Date: Fri, 21 Nov 2025 14:06:48 +0900 Subject: [PATCH] Add complex magnitude and phase helpers --- NumFlat/Mat/Mat.Complex.cs | 142 ++++++++++++++++++ NumFlat/Mat/MatrixExtensions.cs | 66 ++++++++ NumFlat/Vec/Vec.Complex.cs | 103 +++++++++++++ NumFlat/Vec/VectorExtensions.cs | 66 ++++++++ NumFlatTest/MatTests/ComplexTests.cs | 81 ++++++++++ NumFlatTest/MatTests/MatrixExtensionsTests.cs | 75 +++++++++ NumFlatTest/VecTests/ComplexTests.cs | 69 +++++++++ NumFlatTest/VecTests/VectorExtensionsTests.cs | 63 ++++++++ 8 files changed, 665 insertions(+) diff --git a/NumFlat/Mat/Mat.Complex.cs b/NumFlat/Mat/Mat.Complex.cs index ed9612f..a3ae033 100644 --- a/NumFlat/Mat/Mat.Complex.cs +++ b/NumFlat/Mat/Mat.Complex.cs @@ -199,5 +199,147 @@ public static void Imaginary(in Mat x, in Mat destination) od += destination.Stride; } } + + /// + /// Computes the magnitude of each element in the complex matrix. + /// + /// + /// The complex matrix. + /// + /// + /// The destination of the magnitudes. + /// + /// + /// This method does not allocate managed heap memory. + /// + public static void Magnitude(in Mat x, in Mat destination) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + ThrowHelper.ThrowIfEmpty(destination, nameof(destination)); + + if (destination.RowCount != x.RowCount) + { + throw new ArgumentException("'destination.RowCount' must match 'x.RowCount'."); + } + + if (destination.ColCount != x.ColCount) + { + throw new ArgumentException("'destination.ColCount' must match 'x.ColCount'."); + } + + var sx = x.Memory.Span; + var sd = destination.Memory.Span; + var ox = 0; + var od = 0; + while (od < sd.Length) + { + var px = ox; + var pd = od; + var end = od + destination.RowCount; + while (pd < end) + { + sd[pd] = sx[px].Magnitude; + px++; + pd++; + } + ox += x.Stride; + od += destination.Stride; + } + } + + /// + /// Computes the squared magnitude of each element in the complex matrix. + /// + /// + /// The complex matrix. + /// + /// + /// The destination of the squared magnitudes. + /// + /// + /// This method does not allocate managed heap memory. + /// + public static void MagnitudeSquared(in Mat x, in Mat destination) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + ThrowHelper.ThrowIfEmpty(destination, nameof(destination)); + + if (destination.RowCount != x.RowCount) + { + throw new ArgumentException("'destination.RowCount' must match 'x.RowCount'."); + } + + if (destination.ColCount != x.ColCount) + { + throw new ArgumentException("'destination.ColCount' must match 'x.ColCount'."); + } + + var sx = x.Memory.Span; + var sd = destination.Memory.Span; + var ox = 0; + var od = 0; + while (od < sd.Length) + { + var px = ox; + var pd = od; + var end = od + destination.RowCount; + while (pd < end) + { + var value = sx[px]; + sd[pd] = value.Real * value.Real + value.Imaginary * value.Imaginary; + px++; + pd++; + } + ox += x.Stride; + od += destination.Stride; + } + } + + /// + /// Computes the phase of each element in the complex matrix. + /// + /// + /// The complex matrix. + /// + /// + /// The destination of the phases. + /// + /// + /// This method does not allocate managed heap memory. + /// + public static void Phase(in Mat x, in Mat destination) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + ThrowHelper.ThrowIfEmpty(destination, nameof(destination)); + + if (destination.RowCount != x.RowCount) + { + throw new ArgumentException("'destination.RowCount' must match 'x.RowCount'."); + } + + if (destination.ColCount != x.ColCount) + { + throw new ArgumentException("'destination.ColCount' must match 'x.ColCount'."); + } + + var sx = x.Memory.Span; + var sd = destination.Memory.Span; + var ox = 0; + var od = 0; + while (od < sd.Length) + { + var px = ox; + var pd = od; + var end = od + destination.RowCount; + while (pd < end) + { + sd[pd] = sx[px].Phase; + px++; + pd++; + } + ox += x.Stride; + od += destination.Stride; + } + } } } diff --git a/NumFlat/Mat/MatrixExtensions.cs b/NumFlat/Mat/MatrixExtensions.cs index 32519d9..d8e4af0 100644 --- a/NumFlat/Mat/MatrixExtensions.cs +++ b/NumFlat/Mat/MatrixExtensions.cs @@ -212,6 +212,72 @@ public static Mat Imaginary(in this Mat x) return result; } + /// + /// Computes the magnitude of each element in the complex matrix. + /// + /// + /// The complex matrix. + /// + /// + /// The magnitudes of the complex matrix. + /// + /// + /// This method allocates a new matrix which is independent from the original matrix. + /// To avoid the allocation, use instead. + /// + public static Mat Magnitude(in this Mat x) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + + var result = new Mat(x.RowCount, x.ColCount); + Mat.Magnitude(x, result); + return result; + } + + /// + /// Computes the squared magnitude of each element in the complex matrix. + /// + /// + /// The complex matrix. + /// + /// + /// The squared magnitudes of the complex matrix. + /// + /// + /// This method allocates a new matrix which is independent from the original matrix. + /// To avoid the allocation, use instead. + /// + public static Mat MagnitudeSquared(in this Mat x) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + + var result = new Mat(x.RowCount, x.ColCount); + Mat.MagnitudeSquared(x, result); + return result; + } + + /// + /// Computes the phase of each element in the complex matrix. + /// + /// + /// The complex matrix. + /// + /// + /// The phases of the complex matrix. + /// + /// + /// This method allocates a new matrix which is independent from the original matrix. + /// To avoid the allocation, use instead. + /// + public static Mat Phase(in this Mat x) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + + var result = new Mat(x.RowCount, x.ColCount); + Mat.Phase(x, result); + return result; + } + /// /// Computes a matrix inversion, X^-1. /// diff --git a/NumFlat/Vec/Vec.Complex.cs b/NumFlat/Vec/Vec.Complex.cs index c0c2778..5d47c4f 100644 --- a/NumFlat/Vec/Vec.Complex.cs +++ b/NumFlat/Vec/Vec.Complex.cs @@ -102,5 +102,108 @@ public static void Imaginary(in Vec x, in Vec destination) pd += destination.Stride; } } + + /// + /// Computes the magnitude of each element in the complex vector. + /// + /// + /// The complex vector. + /// + /// + /// The destination of the magnitudes. + /// + /// + /// This method does not allocate managed heap memory. + /// + public static void Magnitude(in Vec x, in Vec destination) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + ThrowHelper.ThrowIfEmpty(destination, nameof(destination)); + + if (destination.Count != x.Count) + { + throw new ArgumentException("'destination.Count' must match 'x.Count'."); + } + + var sx = x.Memory.Span; + var sd = destination.Memory.Span; + var px = 0; + var pd = 0; + while (pd < sd.Length) + { + sd[pd] = sx[px].Magnitude; + px += x.Stride; + pd += destination.Stride; + } + } + + /// + /// Computes the squared magnitude of each element in the complex vector. + /// + /// + /// The complex vector. + /// + /// + /// The destination of the squared magnitudes. + /// + /// + /// This method does not allocate managed heap memory. + /// + public static void MagnitudeSquared(in Vec x, in Vec destination) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + ThrowHelper.ThrowIfEmpty(destination, nameof(destination)); + + if (destination.Count != x.Count) + { + throw new ArgumentException("'destination.Count' must match 'x.Count'."); + } + + var sx = x.Memory.Span; + var sd = destination.Memory.Span; + var px = 0; + var pd = 0; + while (pd < sd.Length) + { + var value = sx[px]; + sd[pd] = value.Real * value.Real + value.Imaginary * value.Imaginary; + px += x.Stride; + pd += destination.Stride; + } + } + + /// + /// Computes the phase of each element in the complex vector. + /// + /// + /// The complex vector. + /// + /// + /// The destination of the phases. + /// + /// + /// This method does not allocate managed heap memory. + /// + public static void Phase(in Vec x, in Vec destination) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + ThrowHelper.ThrowIfEmpty(destination, nameof(destination)); + + if (destination.Count != x.Count) + { + throw new ArgumentException("'destination.Count' must match 'x.Count'."); + } + + var sx = x.Memory.Span; + var sd = destination.Memory.Span; + var px = 0; + var pd = 0; + while (pd < sd.Length) + { + sd[pd] = sx[px].Phase; + px += x.Stride; + pd += destination.Stride; + } + } } } diff --git a/NumFlat/Vec/VectorExtensions.cs b/NumFlat/Vec/VectorExtensions.cs index dda2593..2ddf2cf 100644 --- a/NumFlat/Vec/VectorExtensions.cs +++ b/NumFlat/Vec/VectorExtensions.cs @@ -264,6 +264,72 @@ public static Vec Imaginary(in this Vec x) return result; } + /// + /// Computes the magnitude of each element in the complex vector. + /// + /// + /// The complex vector. + /// + /// + /// The magnitudes of the complex vector. + /// + /// + /// This method allocates a new vector which is independent from the original vector. + /// To avoid the allocation, use instead. + /// + public static Vec Magnitude(in this Vec x) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + + var result = new Vec(x.Count); + Vec.Magnitude(x, result); + return result; + } + + /// + /// Computes the squared magnitude of each element in the complex vector. + /// + /// + /// The complex vector. + /// + /// + /// The squared magnitudes of the complex vector. + /// + /// + /// This method allocates a new vector which is independent from the original vector. + /// To avoid the allocation, use instead. + /// + public static Vec MagnitudeSquared(in this Vec x) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + + var result = new Vec(x.Count); + Vec.MagnitudeSquared(x, result); + return result; + } + + /// + /// Computes the phase of each element in the complex vector. + /// + /// + /// The complex vector. + /// + /// + /// The phases of the complex vector. + /// + /// + /// This method allocates a new vector which is independent from the original vector. + /// To avoid the allocation, use instead. + /// + public static Vec Phase(in this Vec x) + { + ThrowHelper.ThrowIfEmpty(x, nameof(x)); + + var result = new Vec(x.Count); + Vec.Phase(x, result); + return result; + } + /// /// Converts the vector to a single-row matrix. /// diff --git a/NumFlatTest/MatTests/ComplexTests.cs b/NumFlatTest/MatTests/ComplexTests.cs index f6bcf00..cca67e3 100644 --- a/NumFlatTest/MatTests/ComplexTests.cs +++ b/NumFlatTest/MatTests/ComplexTests.cs @@ -89,5 +89,86 @@ public void Imaginary(int rowCount, int colCount, int xStride, int dstStride) TestMatrix.FailIfOutOfRangeWrite(dst); } + + [TestCase(1, 1, 1, 1)] + [TestCase(1, 1, 3, 5)] + [TestCase(2, 2, 2, 2)] + [TestCase(2, 2, 4, 3)] + [TestCase(3, 2, 3, 3)] + [TestCase(3, 2, 4, 5)] + [TestCase(2, 3, 2, 3)] + [TestCase(2, 3, 4, 4)] + [TestCase(6, 3, 8, 7)] + [TestCase(4, 7, 5, 9)] + public void Magnitude(int rowCount, int colCount, int xStride, int dstStride) + { + var x = TestMatrix.RandomComplex(42, rowCount, colCount, xStride); + var dst = TestMatrix.RandomDouble(0, rowCount, colCount, dstStride); + + using (x.EnsureUnchanged()) + { + Mat.Magnitude(x, dst); + } + + var expected = x.Map(value => value.Magnitude); + + NumAssert.AreSame(expected, dst, 0); + + TestMatrix.FailIfOutOfRangeWrite(dst); + } + + [TestCase(1, 1, 1, 1)] + [TestCase(1, 1, 3, 5)] + [TestCase(2, 2, 2, 2)] + [TestCase(2, 2, 4, 3)] + [TestCase(3, 2, 3, 3)] + [TestCase(3, 2, 4, 5)] + [TestCase(2, 3, 2, 3)] + [TestCase(2, 3, 4, 4)] + [TestCase(6, 3, 8, 7)] + [TestCase(4, 7, 5, 9)] + public void MagnitudeSquared(int rowCount, int colCount, int xStride, int dstStride) + { + var x = TestMatrix.RandomComplex(42, rowCount, colCount, xStride); + var dst = TestMatrix.RandomDouble(0, rowCount, colCount, dstStride); + + using (x.EnsureUnchanged()) + { + Mat.MagnitudeSquared(x, dst); + } + + var expected = x.Map(value => value.Real * value.Real + value.Imaginary * value.Imaginary); + + NumAssert.AreSame(expected, dst, 0); + + TestMatrix.FailIfOutOfRangeWrite(dst); + } + + [TestCase(1, 1, 1, 1)] + [TestCase(1, 1, 3, 5)] + [TestCase(2, 2, 2, 2)] + [TestCase(2, 2, 4, 3)] + [TestCase(3, 2, 3, 3)] + [TestCase(3, 2, 4, 5)] + [TestCase(2, 3, 2, 3)] + [TestCase(2, 3, 4, 4)] + [TestCase(6, 3, 8, 7)] + [TestCase(4, 7, 5, 9)] + public void Phase(int rowCount, int colCount, int xStride, int dstStride) + { + var x = TestMatrix.RandomComplex(42, rowCount, colCount, xStride); + var dst = TestMatrix.RandomDouble(0, rowCount, colCount, dstStride); + + using (x.EnsureUnchanged()) + { + Mat.Phase(x, dst); + } + + var expected = x.Map(value => value.Phase); + + NumAssert.AreSame(expected, dst, 0); + + TestMatrix.FailIfOutOfRangeWrite(dst); + } } } diff --git a/NumFlatTest/MatTests/MatrixExtensionsTests.cs b/NumFlatTest/MatTests/MatrixExtensionsTests.cs index 47f707e..13e0d57 100644 --- a/NumFlatTest/MatTests/MatrixExtensionsTests.cs +++ b/NumFlatTest/MatTests/MatrixExtensionsTests.cs @@ -231,6 +231,81 @@ public void Imaginary(int rowCount, int colCount, int xStride) NumAssert.AreSame(expected, actual, 0); } + [TestCase(1, 1, 1)] + [TestCase(2, 2, 2)] + [TestCase(2, 2, 4)] + [TestCase(3, 1, 3)] + [TestCase(1, 3, 1)] + [TestCase(1, 3, 2)] + [TestCase(3, 2, 3)] + [TestCase(2, 3, 2)] + [TestCase(4, 5, 8)] + [TestCase(9, 6, 11)] + public void Magnitude(int rowCount, int colCount, int xStride) + { + var x = TestMatrix.RandomComplex(42, rowCount, colCount, xStride); + + Mat actual; + using (x.EnsureUnchanged()) + { + actual = x.Magnitude(); + } + + var expected = x.Map(value => value.Magnitude); + + NumAssert.AreSame(expected, actual, 0); + } + + [TestCase(1, 1, 1)] + [TestCase(2, 2, 2)] + [TestCase(2, 2, 4)] + [TestCase(3, 1, 3)] + [TestCase(1, 3, 1)] + [TestCase(1, 3, 2)] + [TestCase(3, 2, 3)] + [TestCase(2, 3, 2)] + [TestCase(4, 5, 8)] + [TestCase(9, 6, 11)] + public void MagnitudeSquared(int rowCount, int colCount, int xStride) + { + var x = TestMatrix.RandomComplex(42, rowCount, colCount, xStride); + + Mat actual; + using (x.EnsureUnchanged()) + { + actual = x.MagnitudeSquared(); + } + + var expected = x.Map(value => value.Real * value.Real + value.Imaginary * value.Imaginary); + + NumAssert.AreSame(expected, actual, 0); + } + + [TestCase(1, 1, 1)] + [TestCase(2, 2, 2)] + [TestCase(2, 2, 4)] + [TestCase(3, 1, 3)] + [TestCase(1, 3, 1)] + [TestCase(1, 3, 2)] + [TestCase(3, 2, 3)] + [TestCase(2, 3, 2)] + [TestCase(4, 5, 8)] + [TestCase(9, 6, 11)] + public void Phase(int rowCount, int colCount, int xStride) + { + var x = TestMatrix.RandomComplex(42, rowCount, colCount, xStride); + + Mat actual; + using (x.EnsureUnchanged()) + { + actual = x.Phase(); + } + + var expected = x.Map(value => value.Phase); + + NumAssert.AreSame(expected, actual, 0); + } + [TestCase(1, 1)] [TestCase(2, 2)] [TestCase(3, 3)] diff --git a/NumFlatTest/VecTests/ComplexTests.cs b/NumFlatTest/VecTests/ComplexTests.cs index 6ebe417..af24884 100644 --- a/NumFlatTest/VecTests/ComplexTests.cs +++ b/NumFlatTest/VecTests/ComplexTests.cs @@ -77,5 +77,74 @@ public void Imaginary(int count, int xStride, int dstStride) TestVector.FailIfOutOfRangeWrite(actual); } + + [TestCase(1, 1, 1)] + [TestCase(2, 2, 2)] + [TestCase(3, 3, 3)] + [TestCase(1, 3, 4)] + [TestCase(2, 5, 4)] + [TestCase(5, 7, 6)] + public void Magnitude(int count, int xStride, int dstStride) + { + var x = TestVector.RandomComplex(42, count, xStride); + + var expected = x.Select(c => c.Magnitude).ToVector(); + + var actual = TestVector.RandomDouble(0, count, dstStride); + using (x.EnsureUnchanged()) + { + Vec.Magnitude(x, actual); + } + + NumAssert.AreSame(expected, actual, 0); + + TestVector.FailIfOutOfRangeWrite(actual); + } + + [TestCase(1, 1, 1)] + [TestCase(2, 2, 2)] + [TestCase(3, 3, 3)] + [TestCase(1, 3, 4)] + [TestCase(2, 5, 4)] + [TestCase(5, 7, 6)] + public void MagnitudeSquared(int count, int xStride, int dstStride) + { + var x = TestVector.RandomComplex(42, count, xStride); + + var expected = x.Select(c => c.Real * c.Real + c.Imaginary * c.Imaginary).ToVector(); + + var actual = TestVector.RandomDouble(0, count, dstStride); + using (x.EnsureUnchanged()) + { + Vec.MagnitudeSquared(x, actual); + } + + NumAssert.AreSame(expected, actual, 0); + + TestVector.FailIfOutOfRangeWrite(actual); + } + + [TestCase(1, 1, 1)] + [TestCase(2, 2, 2)] + [TestCase(3, 3, 3)] + [TestCase(1, 3, 4)] + [TestCase(2, 5, 4)] + [TestCase(5, 7, 6)] + public void Phase(int count, int xStride, int dstStride) + { + var x = TestVector.RandomComplex(42, count, xStride); + + var expected = x.Select(c => c.Phase).ToVector(); + + var actual = TestVector.RandomDouble(0, count, dstStride); + using (x.EnsureUnchanged()) + { + Vec.Phase(x, actual); + } + + NumAssert.AreSame(expected, actual, 0); + + TestVector.FailIfOutOfRangeWrite(actual); + } } } diff --git a/NumFlatTest/VecTests/VectorExtensionsTests.cs b/NumFlatTest/VecTests/VectorExtensionsTests.cs index ee3d806..5be0e0e 100644 --- a/NumFlatTest/VecTests/VectorExtensionsTests.cs +++ b/NumFlatTest/VecTests/VectorExtensionsTests.cs @@ -278,6 +278,69 @@ public void Imaginary(int count, int xStride) NumAssert.AreSame(expected, actual, 0); } + [TestCase(1, 1)] + [TestCase(2, 2)] + [TestCase(3, 3)] + [TestCase(1, 3)] + [TestCase(2, 5)] + [TestCase(5, 7)] + public void Magnitude(int count, int xStride) + { + var x = TestVector.RandomComplex(42, count, xStride); + + Vec actual; + using (x.EnsureUnchanged()) + { + actual = x.Magnitude(); + } + + var expected = x.Select(value => value.Magnitude).ToVector(); + + NumAssert.AreSame(expected, actual, 0); + } + + [TestCase(1, 1)] + [TestCase(2, 2)] + [TestCase(3, 3)] + [TestCase(1, 3)] + [TestCase(2, 5)] + [TestCase(5, 7)] + public void MagnitudeSquared(int count, int xStride) + { + var x = TestVector.RandomComplex(42, count, xStride); + + Vec actual; + using (x.EnsureUnchanged()) + { + actual = x.MagnitudeSquared(); + } + + var expected = x.Select(value => value.Real * value.Real + value.Imaginary * value.Imaginary).ToVector(); + + NumAssert.AreSame(expected, actual, 0); + } + + [TestCase(1, 1)] + [TestCase(2, 2)] + [TestCase(3, 3)] + [TestCase(1, 3)] + [TestCase(2, 5)] + [TestCase(5, 7)] + public void Phase(int count, int xStride) + { + var x = TestVector.RandomComplex(42, count, xStride); + + Vec actual; + using (x.EnsureUnchanged()) + { + actual = x.Phase(); + } + + var expected = x.Select(value => value.Phase).ToVector(); + + NumAssert.AreSame(expected, actual, 0); + } + [TestCase(1, 1)] [TestCase(2, 2)] [TestCase(3, 3)]