tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Merge pull request #366 from MaLiN2223/master 

More accurate algorithms for Complex
pull/367/head
Christoph Ruegg 11 years ago
parent
3e15a2ac89 6d16f809e0
commit
aad5c5303e
2 changed files with 130 additions and 14 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 85
      src/Numerics/Complex32.cs
  2. 59
      src/UnitTests/ComplexTests/Complex32Test.cs

85
src/Numerics/Complex32.cs
View File

@ -180,11 +180,35 @@ namespace MathNet.Numerics
/// <summary>
/// Gets the magnitude (or absolute value) of a complex number.
/// </summary>
/// <remarks>Assuming that magnitude of (inf,a) and (a,inf) and (inf,inf) is inf and (NaN,a), (a,NaN) and (NaN,NaN) is NaN</remarks>
/// <returns>The magnitude of the current instance.</returns>
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
{
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 (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);
}
}
}
/// <summary>
@ -485,9 +509,9 @@ namespace MathNet.Numerics
/// </summary>
public Tuple<Complex32, Complex32, Complex32> 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<Complex32, Complex32, Complex32>(
FromPolarCoordinates(r, theta),
FromPolarCoordinates(r, theta + shift),
@ -620,6 +644,8 @@ namespace MathNet.Numerics
}
/// <summary>Division operator. Divides a complex number by another.</summary>
/// <remarks>Enhanced Smith's algorithm for dividing two complex numbers </remarks>
/// <see cref="InternalDiv(float, float, float, float, bool)"/>
/// <returns>The result of the division.</returns>
/// <param name="dividend">The dividend.</param>
/// <param name="divisor">The divisor.</param>
@ -634,14 +660,46 @@ 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);
}
/// <summary>
/// Helper method for dividing.
/// </summary>
/// <param name="a">Re first</param>
/// <param name="b">Im first</param>
/// <param name="c">Re second</param>
/// <param name="d">Im second</param>
/// <param name="swapped"></param>
/// <returns></returns>
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);
}
/// <summary>Division operator. Divides a float value by a complex number.</summary>
/// <remarks>Algorithm based on Smith's algorithm</remarks>
/// <see cref="InternalDiv(float, float, float, float, bool)"/>
/// <returns>The result of the division.</returns>
/// <param name="dividend">The dividend.</param>
/// <param name="divisor">The divisor.</param>
@ -656,9 +714,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);
}
/// <summary>Division operator. Divides a complex number by a float value.</summary>
@ -1272,7 +1332,6 @@ namespace MathNet.Numerics
{
return dividend / divisor;
}
/// <summary>
/// Returns the multiplicative inverse of a complex number.
/// </summary>

59
src/UnitTests/ComplexTests/Complex32Test.cs
View File

@ -31,7 +31,7 @@ namespace MathNet.Numerics.UnitTests.ComplexTests
{
#if NOSYSNUMERICS
using Complex = Numerics.Complex;
#else
#else
using Complex = System.Numerics.Complex;
#endif
@ -417,7 +417,54 @@ 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);
}
/// <summary>
/// Can divide without overflow.
/// </summary>
[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);
}
/// <summary>
/// Can divide float/complex without overflow
/// </summary>
[Test]
public void CanDodgeOverflowDivisionFloat()
{
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);
}
/// <summary>
/// Can multiple a complex number and a double using operators.
/// </summary>
@ -547,7 +594,17 @@ namespace MathNet.Numerics.UnitTests.ComplexTests
{
Assert.AreEqual(expected, new Complex32(real, imag).Magnitude);
}
/// <summary>
/// Can calculate magnitude without overflow
/// </summary>
[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);
}
/// <summary>
/// Can compute sign.
/// </summary>

Write Preview
Loading…
Cancel
Save