Changed dividing and magintude algorithms for more accurate.

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);
+                } 
+            }
         }
 
         /// <summary>
@@ -485,9 +504,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 +639,7 @@ namespace MathNet.Numerics
         }
 
         /// <summary>Division operator. Divides a complex number by another.</summary>
+        /// <remarks>Enchanted Smith's algorithm for dividing two complex numbers </remarks>
         /// <returns>The result of the division.</returns>
         /// <param name="dividend">The dividend.</param>
         /// <param name="divisor">The divisor.</param>
@@ -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);
+        }
+        /// <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>
@@ -1272,7 +1322,6 @@ namespace MathNet.Numerics
         {
             return dividend / divisor;
         }
-
         /// <summary>
         /// Returns the multiplicative inverse of a complex number.
         /// </summary>

src/UnitTests/ComplexTests/Complex32Test.cs

@@ -24,14 +24,15 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
-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);
         }
+        /// <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 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);
         }
+        /// <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>