diff --git a/src/Numerics/Complex32.cs b/src/Numerics/Complex32.cs
index a134cc38..519bab95 100644
--- a/src/Numerics/Complex32.cs
+++ b/src/Numerics/Complex32.cs
@@ -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>
diff --git a/src/UnitTests/ComplexTests/Complex32Test.cs b/src/UnitTests/ComplexTests/Complex32Test.cs
index ec28709c..c7f20678 100644
--- a/src/UnitTests/ComplexTests/Complex32Test.cs
+++ b/src/UnitTests/ComplexTests/Complex32Test.cs
@@ -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>