diff --git a/src/Numerics/ComplexExtensions.cs b/src/Numerics/ComplexExtensions.cs index 4d8f8487..8050e43d 100644 --- a/src/Numerics/ComplexExtensions.cs +++ b/src/Numerics/ComplexExtensions.cs @@ -243,10 +243,47 @@ namespace MathNet.Numerics /// /// The square root of this complex number. /// - [TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")] public static Complex SquareRoot(this Complex complex) { - return Complex.Sqrt(complex); + // Note: the following code should be equivalent to Complex.Sqrt(complex), + // but it turns out that is implemented poorly in System.Numerics, + // hence we provide our own implementation here. Do not replace. + + if (complex.IsRealNonNegative()) + { + return new Complex(Math.Sqrt(complex.Real), 0.0); + } + + Complex result; + + var absReal = Math.Abs(complex.Real); + var absImag = Math.Abs(complex.Imaginary); + double w; + if (absReal >= absImag) + { + var ratio = complex.Imaginary / complex.Real; + w = Math.Sqrt(absReal) * Math.Sqrt(0.5 * (1.0 + Math.Sqrt(1.0 + (ratio * ratio)))); + } + else + { + var ratio = complex.Real / complex.Imaginary; + w = Math.Sqrt(absImag) * Math.Sqrt(0.5 * (Math.Abs(ratio) + Math.Sqrt(1.0 + (ratio * ratio)))); + } + + if (complex.Real >= 0.0) + { + result = new Complex(w, complex.Imaginary / (2.0 * w)); + } + else if (complex.Imaginary >= 0.0) + { + result = new Complex(absImag / (2.0 * w), w); + } + else + { + result = new Complex(absImag / (2.0 * w), -w); + } + + return result; } ///