From 5b9aced350e359e6c4ccede898159b6ab8ecf84c Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Sat, 23 Jan 2016 17:33:11 +0100 Subject: [PATCH 1/9] Changed dividing and magintude algorithms for more accurate. --- src/Numerics/Complex32.cs | 69 ++++++++++++++++++--- src/UnitTests/ComplexTests/Complex32Test.cs | 32 +++++++++- 2 files changed, 90 insertions(+), 11 deletions(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index a134cc38..39b7de51 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -184,7 +184,26 @@ namespace MathNet.Numerics public float Magnitude { [TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")] - get { return (float)Math.Sqrt((_real * _real) + (_imag * _imag)); } + get + { + float a = Math.Abs(_real); + float b = Math.Abs(_imag); + if (a > b) + { + double tmp = b / a; + return a * (float)Math.Sqrt(1.0f + tmp * tmp); + + } + if (a == 0.0f) // one can write a >= float.Epsilon here + { + return b; + } + else + { + double tmp = a / b; + return b*(float)Math.Sqrt(1.0f + tmp*tmp); + } + } } /// @@ -485,9 +504,9 @@ namespace MathNet.Numerics /// public Tuple CubicRoots() { - float r = (float)Math.Pow(Magnitude, 1d/3d); - float theta = Phase/3; - const float shift = (float)Constants.Pi2/3; + float r = (float)Math.Pow(Magnitude, 1d / 3d); + float theta = Phase / 3; + const float shift = (float)Constants.Pi2 / 3; return new Tuple( FromPolarCoordinates(r, theta), FromPolarCoordinates(r, theta + shift), @@ -620,6 +639,7 @@ namespace MathNet.Numerics } /// Division operator. Divides a complex number by another. + /// Enchanted Smith's algorithm for dividing two complex numbers /// The result of the division. /// The dividend. /// The divisor. @@ -634,11 +654,41 @@ namespace MathNet.Numerics { return PositiveInfinity; } - - var modSquared = divisor.MagnitudeSquared; - return new Complex32( - ((dividend._real * divisor._real) + (dividend._imag * divisor._imag)) / modSquared, - ((dividend._imag * divisor._real) - (dividend._real * divisor._imag)) / modSquared); + float a = dividend.Real; + float b = dividend.Imaginary; + float c = divisor.Real; + float d = divisor.Imaginary; + if (Math.Abs(d) <= Math.Abs(c)) + return InternalDiv(a, b, c, d, false); + return InternalDiv(b, a, d, c, true); + } + /// + /// Helper method for dividing. + /// + /// Re first + /// Im first + /// Re second + /// Im second + /// + /// + private static Complex32 InternalDiv(float a, float b, float c, float d, bool swapped) + { + float r = d / c; + float t = 1 / (c + d * r); + float e, f; + if (r!=0.0f) // one can use r >= float.Epsilon || r <= float.Epsilon instead + { + e = (a + b * r) * t; + f = (b - a * r) * t; + } + else + { + e = (a + d * (b / c)) * t; + f = (b - d * (a / c)) * t; + } + if (swapped) + f = -f; + return new Complex32(e, f); } /// Division operator. Divides a float value by a complex number. @@ -1272,7 +1322,6 @@ namespace MathNet.Numerics { return dividend / divisor; } - /// /// Returns the multiplicative inverse of a complex number. /// diff --git a/src/UnitTests/ComplexTests/Complex32Test.cs b/src/UnitTests/ComplexTests/Complex32Test.cs index ec28709c..fd5f131b 100644 --- a/src/UnitTests/ComplexTests/Complex32Test.cs +++ b/src/UnitTests/ComplexTests/Complex32Test.cs @@ -24,14 +24,15 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; using NUnit.Framework; +using System; namespace MathNet.Numerics.UnitTests.ComplexTests { #if NOSYSNUMERICS using Complex = Numerics.Complex; #else + using System.Diagnostics; using Complex = System.Numerics.Complex; #endif @@ -417,6 +418,25 @@ namespace MathNet.Numerics.UnitTests.ComplexTests Assert.AreEqual(new Complex32(-2, 0), new Complex32(4, -4) / new Complex32(-2, 2)); Assert.AreEqual(Complex32.PositiveInfinity, Complex32.One / Complex32.Zero); } + /// + /// Can divide without overflow. + /// + [Test] + public void CanDodgeOverflowDivision() + { + var first = new Complex32((float)Math.Pow(10, 37), (float)Math.Pow(10, -37)); + var second = new Complex32((float)Math.Pow(10, 25), (float)Math.Pow(10, -25)); + Assert.AreEqual(new Complex32((float)Math.Pow(10, 12), (float)Math.Pow(10, -38)), first / second); + + first = new Complex32(-(float)Math.Pow(10, 37), (float)Math.Pow(10, -37)); + second = new Complex32((float)Math.Pow(10, 25), (float)Math.Pow(10, -25)); + Assert.AreEqual(new Complex32(-(float)Math.Pow(10, 12), (float)Math.Pow(10, -38)), first / second); + + first = new Complex32((float)Math.Pow(10, -37), (float)Math.Pow(10, 37)); + second = new Complex32((float)Math.Pow(10, -17), -(float)Math.Pow(10, 17)); + Assert.AreEqual(new Complex32(-(float)Math.Pow(10, 20), (float)Math.Pow(10, -14)), first / second); + + } /// /// Can multiple a complex number and a double using operators. @@ -547,7 +567,17 @@ namespace MathNet.Numerics.UnitTests.ComplexTests { Assert.AreEqual(expected, new Complex32(real, imag).Magnitude); } + /// + /// Can calculate magnitude without overflow + /// + [Test] + public void CanDodgeOverflowMagnitude() + { + Assert.AreEqual((float)Math.Sqrt(2) * float.Epsilon, new Complex32(float.Epsilon, float.Epsilon).Magnitude); + Assert.AreEqual(float.Epsilon, new Complex32(0, float.Epsilon).Magnitude); + Assert.AreEqual((float)(Math.Pow(10,30) * Math.Sqrt(2)), new Complex32((float)Math.Pow(10, 30), (float)Math.Pow(10, 30)).Magnitude); + } /// /// Can compute sign. /// From 30132b285cf737b9d5eb00a6bb1903e3db919e03 Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Sat, 23 Jan 2016 17:36:49 +0100 Subject: [PATCH 2/9] formatting changes --- src/Numerics/Complex32.cs | 6 +++--- src/UnitTests/ComplexTests/Complex32Test.cs | 5 ++--- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index 39b7de51..2024bb60 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -504,9 +504,9 @@ namespace MathNet.Numerics /// public Tuple CubicRoots() { - float r = (float)Math.Pow(Magnitude, 1d / 3d); - float theta = Phase / 3; - const float shift = (float)Constants.Pi2 / 3; + float r = (float)Math.Pow(Magnitude, 1d/3d); + float theta = Phase/3; + const float shift = (float)Constants.Pi2/3; return new Tuple( FromPolarCoordinates(r, theta), FromPolarCoordinates(r, theta + shift), diff --git a/src/UnitTests/ComplexTests/Complex32Test.cs b/src/UnitTests/ComplexTests/Complex32Test.cs index fd5f131b..d381a146 100644 --- a/src/UnitTests/ComplexTests/Complex32Test.cs +++ b/src/UnitTests/ComplexTests/Complex32Test.cs @@ -24,15 +24,14 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using NUnit.Framework; using System; +using NUnit.Framework; namespace MathNet.Numerics.UnitTests.ComplexTests { #if NOSYSNUMERICS using Complex = Numerics.Complex; -#else - using System.Diagnostics; +#else using Complex = System.Numerics.Complex; #endif From 9382d19f00c17d8ea0ee06e3f65a0c777ff4a4bf Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Sat, 23 Jan 2016 18:23:50 +0100 Subject: [PATCH 3/9] Magnitude of Complex(inf,inf) == inf --- src/Numerics/Complex32.cs | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index 2024bb60..6ccd0de3 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -188,6 +188,8 @@ namespace MathNet.Numerics { float a = Math.Abs(_real); float b = Math.Abs(_imag); + if (float.IsPositiveInfinity(a) && float.IsPositiveInfinity(b)) + return float.PositiveInfinity; if (a > b) { double tmp = b / a; @@ -201,8 +203,8 @@ namespace MathNet.Numerics else { double tmp = a / b; - return b*(float)Math.Sqrt(1.0f + tmp*tmp); - } + return b * (float)Math.Sqrt(1.0f + tmp * tmp); + } } } @@ -504,9 +506,9 @@ namespace MathNet.Numerics /// public Tuple CubicRoots() { - float r = (float)Math.Pow(Magnitude, 1d/3d); - float theta = Phase/3; - const float shift = (float)Constants.Pi2/3; + float r = (float)Math.Pow(Magnitude, 1d / 3d); + float theta = Phase / 3; + const float shift = (float)Constants.Pi2 / 3; return new Tuple( FromPolarCoordinates(r, theta), FromPolarCoordinates(r, theta + shift), @@ -676,7 +678,7 @@ namespace MathNet.Numerics float r = d / c; float t = 1 / (c + d * r); float e, f; - if (r!=0.0f) // one can use r >= float.Epsilon || r <= float.Epsilon instead + if (r != 0.0f) // one can use r >= float.Epsilon || r <= float.Epsilon instead { e = (a + b * r) * t; f = (b - a * r) * t; From 769764464c62d9ef1f16768b38282e88502efa48 Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Sun, 24 Jan 2016 13:41:41 +0100 Subject: [PATCH 4/9] Corrected typo, corrected float/complex, added more tests. --- src/Numerics/Complex32.cs | 13 +++++++++---- src/UnitTests/ComplexTests/Complex32Test.cs | 16 +++++++++++++++- 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index 6ccd0de3..18097c2f 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -641,7 +641,8 @@ namespace MathNet.Numerics } /// Division operator. Divides a complex number by another. - /// Enchanted Smith's algorithm for dividing two complex numbers + /// Enhanced Smith's algorithm for dividing two complex numbers + /// /// The result of the division. /// The dividend. /// The divisor. @@ -694,6 +695,8 @@ namespace MathNet.Numerics } /// Division operator. Divides a float value by a complex number. + /// Algorithm based on Smith's algorithm + /// /// The result of the division. /// The dividend. /// The divisor. @@ -708,9 +711,11 @@ namespace MathNet.Numerics { return PositiveInfinity; } - - var zmod = divisor.MagnitudeSquared; - return new Complex32(dividend * divisor._real / zmod, -dividend * divisor._imag / zmod); + float c = divisor.Real; + float d = divisor.Imaginary; + if (Math.Abs(d) <= Math.Abs(c)) + return InternalDiv(dividend, 0, c, d, false); + return InternalDiv(0, dividend, d, c, true); } /// Division operator. Divides a complex number by a float value. diff --git a/src/UnitTests/ComplexTests/Complex32Test.cs b/src/UnitTests/ComplexTests/Complex32Test.cs index d381a146..355f4fe2 100644 --- a/src/UnitTests/ComplexTests/Complex32Test.cs +++ b/src/UnitTests/ComplexTests/Complex32Test.cs @@ -436,7 +436,21 @@ namespace MathNet.Numerics.UnitTests.ComplexTests Assert.AreEqual(new Complex32(-(float)Math.Pow(10, 20), (float)Math.Pow(10, -14)), first / second); } + /// + /// Can divide float/complex without overflow + /// + [Test] + public void CanDodgeOverflowDivisionFloat() + { + + var first = new Complex32((float)Math.Pow(10, 25), (float)Math.Pow(10, -25)); + float second = (float)Math.Pow(10, -37); + float third = (float)Math.Pow(10, 37); + Assert.AreEqual(new Complex32(0, 0), second / first); // it's (10^-62,10^-112) thus overflow + Assert.AreEqual(new Complex32((float)Math.Pow(10, 12), (float)Math.Pow(10, -38)), third / first); + Assert.AreEqual(new Complex32((float)Math.Pow(10, -28), -(float)Math.Pow(10, 12)), third / new Complex32((float)Math.Pow(10, -25), (float)Math.Pow(10, 25))); + } /// /// Can multiple a complex number and a double using operators. /// @@ -574,7 +588,7 @@ namespace MathNet.Numerics.UnitTests.ComplexTests { Assert.AreEqual((float)Math.Sqrt(2) * float.Epsilon, new Complex32(float.Epsilon, float.Epsilon).Magnitude); Assert.AreEqual(float.Epsilon, new Complex32(0, float.Epsilon).Magnitude); - Assert.AreEqual((float)(Math.Pow(10,30) * Math.Sqrt(2)), new Complex32((float)Math.Pow(10, 30), (float)Math.Pow(10, 30)).Magnitude); + Assert.AreEqual((float)(Math.Pow(10, 30) * Math.Sqrt(2)), new Complex32((float)Math.Pow(10, 30), (float)Math.Pow(10, 30)).Magnitude); } /// From 1d71dc53c9a24ecff38864453f50d61d7472363c Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Sun, 24 Jan 2016 13:56:36 +0100 Subject: [PATCH 5/9] || instead of && in complex division. More tests. --- src/Numerics/Complex32.cs | 2 +- src/UnitTests/ComplexTests/Complex32Test.cs | 26 ++++++++++++++++----- 2 files changed, 21 insertions(+), 7 deletions(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index 18097c2f..15e9581d 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -188,7 +188,7 @@ namespace MathNet.Numerics { float a = Math.Abs(_real); float b = Math.Abs(_imag); - if (float.IsPositiveInfinity(a) && float.IsPositiveInfinity(b)) + if (float.IsPositiveInfinity(a) || float.IsPositiveInfinity(b)) return float.PositiveInfinity; if (a > b) { diff --git a/src/UnitTests/ComplexTests/Complex32Test.cs b/src/UnitTests/ComplexTests/Complex32Test.cs index 355f4fe2..c7f20678 100644 --- a/src/UnitTests/ComplexTests/Complex32Test.cs +++ b/src/UnitTests/ComplexTests/Complex32Test.cs @@ -443,12 +443,26 @@ namespace MathNet.Numerics.UnitTests.ComplexTests public void CanDodgeOverflowDivisionFloat() { - var first = new Complex32((float)Math.Pow(10, 25), (float)Math.Pow(10, -25)); - float second = (float)Math.Pow(10, -37); - float third = (float)Math.Pow(10, 37); - Assert.AreEqual(new Complex32(0, 0), second / first); // it's (10^-62,10^-112) thus overflow - Assert.AreEqual(new Complex32((float)Math.Pow(10, 12), (float)Math.Pow(10, -38)), third / first); - Assert.AreEqual(new Complex32((float)Math.Pow(10, -28), -(float)Math.Pow(10, 12)), third / new Complex32((float)Math.Pow(10, -25), (float)Math.Pow(10, 25))); + var firstComplex = new Complex32((float)Math.Pow(10, 25), (float)Math.Pow(10, -25)); + float firstFloat = (float)Math.Pow(10, -37); + float secondFloat = (float)Math.Pow(10, 37); + + Assert.AreEqual(new Complex32(0, 0), firstFloat / firstComplex); // it's (10^-62,10^-112) thus overflow to 0 + Assert.AreEqual(new Complex32((float)Math.Pow(10, 12), (float)Math.Pow(10, -38)), secondFloat / firstComplex); + + var secondComplex = new Complex32((float)Math.Pow(10, -25), (float)Math.Pow(10, 25)); + Assert.AreEqual(new Complex32(0, 0), firstFloat / secondComplex);// it's (10^-112,10^-62) thus overflow to 0 + Assert.AreEqual(new Complex32((float)Math.Pow(10, -38), -(float)Math.Pow(10, 12)), secondFloat / secondComplex); + + + + float thirdFloat = (float)Math.Pow(10, 13); + var thirdComplex = new Complex32((float)Math.Pow(10, -25), (float)Math.Pow(10, -25)); + Assert.AreEqual(new Complex32(5.0f * (float)Math.Pow(10, 37), -5.0f * (float)Math.Pow(10, 37)), thirdFloat / thirdComplex); + + var fourthFloat = (float)Math.Pow(10, -30); + var fourthComplex = new Complex32((float)Math.Pow(10, -30), (float)Math.Pow(10, -30)); + Assert.AreEqual(new Complex32(0.5f, -0.5f), fourthFloat / fourthComplex); } /// From 0f49428b8b8f1d59c9863eb06e6165bbfc79b6fb Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Sun, 24 Jan 2016 14:05:25 +0100 Subject: [PATCH 6/9] Changed || back to && because of MinimumMagnitudePhaseOfNaNIsNaN test in StatisticsTests. --- src/Numerics/Complex32.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index 15e9581d..d2cfd930 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -188,7 +188,7 @@ namespace MathNet.Numerics { float a = Math.Abs(_real); float b = Math.Abs(_imag); - if (float.IsPositiveInfinity(a) || float.IsPositiveInfinity(b)) + if (float.IsPositiveInfinity(a) && float.IsPositiveInfinity(b)) return float.PositiveInfinity; if (a > b) { From 75ba3f82791806383a4864224f3de0b08348f958 Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Sun, 24 Jan 2016 14:07:30 +0100 Subject: [PATCH 7/9] Added remark for Magnitude --- src/Numerics/Complex32.cs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index d2cfd930..608014d0 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -180,6 +180,7 @@ namespace MathNet.Numerics /// /// Gets the magnitude (or absolute value) of a complex number. /// + /// Assuming that magnitude of (inf,a) and (a,inf) is NaN and magnitude of (inf,inf) is inf /// The magnitude of the current instance. public float Magnitude { @@ -188,7 +189,7 @@ namespace MathNet.Numerics { float a = Math.Abs(_real); float b = Math.Abs(_imag); - if (float.IsPositiveInfinity(a) && float.IsPositiveInfinity(b)) + if (float.IsPositiveInfinity(a) && float.IsPositiveInfinity(b)) return float.PositiveInfinity; if (a > b) { From 987704134452650f77c06446dcaec0b4e342b5a7 Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Mon, 25 Jan 2016 11:16:57 +0100 Subject: [PATCH 8/9] added if's in magnitude calculation --- src/Numerics/Complex32.cs | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index 608014d0..cfa0f009 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -187,10 +187,12 @@ namespace MathNet.Numerics [TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")] get { + if (float.IsNaN(_real) || float.IsNaN(_imag)) + return float.NaN; + if (float.IsInfinity(_real) || float.IsInfinity(_imag)) + return float.PositiveInfinity; float a = Math.Abs(_real); float b = Math.Abs(_imag); - if (float.IsPositiveInfinity(a) && float.IsPositiveInfinity(b)) - return float.PositiveInfinity; if (a > b) { double tmp = b / a; From 6d16f809e06f68f3382f2c14a7e735d497c7c5c2 Mon Sep 17 00:00:00 2001 From: MaLiN2223 Date: Mon, 25 Jan 2016 11:49:52 +0100 Subject: [PATCH 9/9] Corrected remark --- src/Numerics/Complex32.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs index cfa0f009..519bab95 100644 --- a/src/Numerics/Complex32.cs +++ b/src/Numerics/Complex32.cs @@ -180,7 +180,7 @@ namespace MathNet.Numerics /// /// Gets the magnitude (or absolute value) of a complex number. /// - /// Assuming that magnitude of (inf,a) and (a,inf) is NaN and magnitude of (inf,inf) is inf + /// Assuming that magnitude of (inf,a) and (a,inf) and (inf,inf) is inf and (NaN,a), (a,NaN) and (NaN,NaN) is NaN /// The magnitude of the current instance. public float Magnitude {