public static Complex[,] GetFilter(double[,] arr)
        {
            int xLength = arr.GetLength(0) - 2 * arr.GetLength(0) / 2 == 0 ? arr.GetLength(0) + 1 : arr.GetLength(0);
            int yLength = arr.GetLength(1) - 2 * arr.GetLength(1) / 2 == 0 ? arr.GetLength(1) + 1 : arr.GetLength(1);
            int upperBoundX = (int)(xLength / 2);
            int lowerBoundX = -1 * upperBoundX;
            int upperBoundY = (int)(yLength / 2);
            int lowerBoundY = -1 * upperBoundY;

            Complex[,] filter = new Complex[xLength, yLength];
            double x_, y_ = 0;

            for (int x = lowerBoundX; x < upperBoundX; x++)
            {
                for (int y = lowerBoundY; y < upperBoundY; y++)
                {
                    if (y < 0)
                    {
                        filter[upperBoundX + x, upperBoundY + y] = 0;
                        continue;
                    }

                    x_ = XY_WithLine.GetX_WithLine(upperBoundX + x, upperBoundY + y);
                    y_ = XY_WithLine.GetY_WithLine(upperBoundX + x, upperBoundY + y);
                    filter[upperBoundX + x, upperBoundY + y] = 
                        (2 * x_ * y_ + Complex.ImaginaryOne * (x_ * x_ - y_ * y_)) * Gaussian.CalculateGaussian(x_, y_);
                }
            }

            return filter;
        }

        public void ForwardInplaceRealSine()
        {
            var samples = Generate.PeriodicMap(16, w => new Complex(Math.Sin(w), 0), 16, 1.0, Constants.Pi2);
            var spectrum = new Complex[samples.Length];

            // real-odd transforms to imaginary odd
            samples.Copy(spectrum);
            Control.FourierTransformProvider.Forward(spectrum, FourierTransformScaling.BackwardScaling);

            // all real components must be zero
            foreach (var c in spectrum)
            {
                Assert.AreEqual(0, c.Real, 1e-12, "real");
            }

            // all imaginary components except second and last musth be zero
            for (var i = 0; i < spectrum.Length; i++)
            {
                if (i == 1)
                {
                    Assert.AreEqual(-8, spectrum[i].Imaginary, 1e-12, "imag second");
                }
                else if (i == spectrum.Length - 1)
                {
                    Assert.AreEqual(8, spectrum[i].Imaginary, 1e-12, "imag last");
                }
                else
                {
                    Assert.AreEqual(0, spectrum[i].Imaginary, 1e-12, "imag");
                }
            }
        }

        /// <summary>
        /// The fourier transform.
        /// </summary>
        /// <param name="characteristics">
        /// The characteristics.
        /// </param>
        /// <returns>
        /// Spectrum of the signal as <see cref="T:double[][]"/>.
        /// </returns>
        public static double[][] CalculateFastFourierTransform(double[][] characteristics)
        {
            var powerOfTwo = PowerOfTwoCeiling(characteristics.Length);
            var result = new double[powerOfTwo][];
            for (int j = 0; j < powerOfTwo; j++)
            {
                result[j] = new double[characteristics[0].Length];
            }

            // transforming into complex representation
            // cycle through all characteristics
            for (int i = 0; i < characteristics[0].Length; i++)
            {
                var complex = new Complex[powerOfTwo];

                // cycle through all sequence fragments
                for (int j = 0; j < powerOfTwo; j++)
                {
                    complex[j] = j < characteristics.Length ? new Complex(characteristics[j][i], 0) : new Complex(0, 0);
                }

                Complex[] complexResult = CalculateFastFourierTransform(complex);

                // converting array to double
                for (int g = 0; g < powerOfTwo; g++)
                {
                    result[g][i] = complexResult[g].Real;
                }
            }

            return result;
        }

        /// <summary>
        /// Recursive solution to determine whether complex number belongs to the Mandelbrot set
        /// </summary>
        public MandelbrotResult IsMemberRecursive(Complex c, Complex? z = null, int iteration = 0)
        {
            if (iteration == MaxIterations) return new MandelbrotResult(iteration, true);
            if (z.HasValue && z.Value.Magnitude > 2.0) return new MandelbrotResult(iteration, false);

            var newZ = z.HasValue ? Complex.Pow(z.Value, 2) + c : Complex.Zero;
            return IsMemberRecursive(c, newZ, iteration + 1);
        }

 public Mobius( Complex a, Complex b, Complex c, Complex d )
     : this()
 {
     A = a;
     B = b;
     C = c;
     D = d;
 }

        public void ComparePowIntWithPow()
        {
            Complex z1 = new Complex(1.6859, 0.3902);
            Complex actual = z1.Pow(10);
            Complex expected = z1 * z1 * z1 * z1 * z1 * z1 * z1 * z1 * z1 * z1;

            Assert.IsTrue(expected.ApproximatelyEquals(actual, epsilon));
        }

        public static PyComplexObject Create(DkmProcess process, Complex value) {
            var allocator = process.GetDataItem<PyObjectAllocator>();
            Debug.Assert(allocator != null);

            var result = allocator.Allocate<PyComplexObject>();
            result.cval.real.Write(value.Real);
            result.cval.imag.Write(value.Imaginary);
            return result;
        }

            public static SparseMatrix PorousMatrix(bool Rigid, double d, Complex k, Complex sin_theta, double freq, double porosity, double tortuosity, double YoungsModulus, double PoissonRatio, double Viscous_Characteristic_Length, double flow_resistivity, double FrameDensity, double Thermal_Permeability_0, double AmbientMeanPressure)
            {
                double w = Utilities.Numerics.PiX2 * freq;
                double v = Biot_Porous_Absorbers.v();
                double FrameShear = AbsorptionModels.Biot_Porous_Absorbers.Shear_Modulus(YoungsModulus, PoissonRatio);
                double kb = 2 * FrameShear * (PoissonRatio + 1) / (3 * (1 - 2 * PoissonRatio));
                double BulkMod_Frame = AbsorptionModels.Biot_Porous_Absorbers.BulkMod_Solid(YoungsModulus, PoissonRatio);
                Complex Kf = Biot_Porous_Absorbers.BulkMod_Fluid(w, AmbientMeanPressure, porosity, Thermal_Permeability_0);//AmbientMeanPressure / (1 - (gamma - 1) / (gamma * alpha));
                Complex LameL = YoungsModulus * PoissonRatio / ((1 + PoissonRatio) * (1 - 2 * PoissonRatio));
                Complex LameMu = YoungsModulus / (2 * (1 + PoissonRatio));
                Complex delta21 = w * w * FrameDensity;
                Complex delta22 = w * w * FrameDensity;
                Complex delta23 = delta21 / LameMu;
                delta21 /= (LameL + 2 * LameMu);
                delta22 /= (LameL + 2 * LameMu);

                //Taken from Lauriks, et. al., 1990.
                double rho12 = Biot_Porous_Absorbers.rho12(porosity, tortuosity);
                double rhoa = Biot_Porous_Absorbers.rhoA(rho12);
                double Viscous_Permeability = Biot_Porous_Absorbers.Viscous_Permeability(flow_resistivity);
                Complex Gw = Biot_Porous_Absorbers.G_w(tortuosity, porosity, Viscous_Permeability, Viscous_Characteristic_Length, freq, v);

                //Complex rho12eff = Biot_Porous_Absorbers.rho12eff(rhoa, porosity, flow_resistivity, Gw, freq);
                Complex rho22eff = Biot_Porous_Absorbers.rho22eff(rhoa, porosity, flow_resistivity, Gw, freq);
                Complex rho11eff = Biot_Porous_Absorbers.rho11eff(FrameDensity, rhoa, porosity, flow_resistivity, Gw, freq);

                Complex P, Q, R;

                if (!Rigid)
                {
                    //Universal (Limp) Frame Case:
                    P = ((1 - porosity) * (1 - kb / BulkMod_Frame) * BulkMod_Frame + porosity * BulkMod_Frame * kb / Kf) / (1 - porosity - kb / BulkMod_Frame + porosity * BulkMod_Frame / Kf);
                    Q = (1 - porosity - kb / BulkMod_Frame) * porosity * BulkMod_Frame / (1 - porosity - kb / BulkMod_Frame + porosity * BulkMod_Frame / Kf);
                    R = porosity * porosity * BulkMod_Frame / (1 - porosity - kb / BulkMod_Frame + porosity * BulkMod_Frame / Kf);
                }
                else
                {
                    //Rigid Frame Case:
                    P = 4 * FrameShear / 3 + kb + (porosity * porosity) * Kf / porosity;
                    R = porosity * Kf;
                    Q = Kf * (1 - porosity);
                }

                Complex kt = k * sin_theta;
                Complex k13 = Complex.Sqrt(delta21 - kt * kt);
                Complex k23 = Complex.Sqrt(delta22 - kt * kt);
                Complex k33 = Complex.Sqrt(delta23 - kt * kt);
                Complex Mu1 = Q * delta21 - w * w * rho11eff / (w * w * rho22eff - R * delta21);
                Complex Mu2 = Q * delta22 - w * w * rho11eff / (w * w * rho22eff - R * delta22);
                Complex Mu3 = FrameShear * delta23 - w * w * rho11eff / (w * w * rho22eff);

                SparseMatrix GH = GammaH_P(kt, w, d, FrameShear, P, Q, R, k13, k23, k33, Mu1, Mu2, Mu3);
                SparseMatrix G0T = Gamma0T_P(kt, w, FrameShear, P, Q, R, k13, k23, k33, Mu1, Mu2, Mu3);

                return GH * G0T;

            }

        public void TestIfMulIsEqualToComplexOperator()
        {
            var z1 = new Complex(1.6859, 0.3902);
            var z2 = new Complex(3.51896, -0.458);

            var w1 = z1 * z2;
            var w2 = ComplexArithmetic.Multiply(z1, z2);

            Assert.IsTrue(w1.ApproximatelyEquals(w2, Epsilon));
        }

        public void TestIfPowIntInFasterThanPow()
        {
            int iter = 250000;
            Complex z = new Complex(1.6859, 0.3902);
            double averageIntTime = BenchmarkUtil.Benchmark(
                () => { Complex w = ComplexArithmetic.PowInt(z, 10); }, iter);
            double averageCpTime = BenchmarkUtil.Benchmark(
                () => { Complex w = Complex.Pow(z, 10); }, iter);

            Assert.IsTrue(averageIntTime < averageCpTime);
        }

 public MandelbrotResult IsMemberIterative(Complex c)
 {
     var iteration = 0;
     var z = Complex.Zero;
     while (iteration < MaxIterations && z.Magnitude < 2)
     {
         z = z * z + c;
         iteration++;
     }
     var isMember = (iteration == MaxIterations);
     return new MandelbrotResult(iteration, isMember);
 }

        public void FourierRadix2IsReversible(FourierOptions options)
        {
            var samples = Generate.RandomComplex(0x8000, GetUniform(1));
            var work = new Complex[samples.Length];
            samples.CopyTo(work, 0);

            Fourier.Radix2Forward(work, options);
            Assert.IsFalse(work.ListAlmostEqual(samples, 6));

            Fourier.Radix2Inverse(work, options);
            AssertHelpers.AlmostEqual(samples, work, 12);
        }

        /// <summary>
        /// Тихоновская регуляризация
        /// </summary>
        /// <param name="filter"> ядро искажения (PSF)</param>
        /// <returns></returns>
        public static ConvolutionFilter Filtering(ConvolutionFilter filter)
        {
            ///в частотной области
            ///fn(u,v)=((h*(u,v)/|h(u,v)|^2+gamma*|p(u,v)|^2))*g(u,v)
            ///fn - приближение
            ///h - kernel
            ///h* - комплексно-сопряженная форма kernel
            ///|h|^2 = h(u,v)*h*(u,v) = u^2+v^2*i
            ///gamma - какой-то параметр (в инверсном фильтре = 0)
            ///p(u,v) = оператор Лапласа = [{0  1  0}
            ///                             {1 -4  1}
            ///                             {0  1  0}]
            ///g - искаженное изображение

            Complex[,] otf = OpticalTransferFunction.Psf2otf(filter);
            int height = otf.GetLength(0);                                              //строк
            int width = otf.GetLength(1);                                              //столбцов
            Complex gamma = Complex.Zero;                                        //
            Complex[,] otfZ = new Complex[height, width];                                   //комплексно сопряженная матрица ядра
            Complex[,] otf2 = new Complex[height, width];                                   //матрица = |h|^2
            Complex[,] p = {{0, 1, 0,},                                          //лапласиан
                           {1, -4, 1,},
                           {0, 1, 0,},};
            p = Fourier.Transform(p);
            for (int u = 0; u < p.GetLength(0); u++)
                for (int v = 0; v < p.GetLength(1); v++)
                    p[u, v] = OpticalTransferFunction.ModPow(p[u, v]);

            for (int u = 0; u < height; u++)
                for (int v = 0; v < width; v++)
                    otfZ[u, v] = Complex.Conjugate(otf[u, v]);

            for (int u = 0; u < height; u++)
                for (int v = 0; v < width; v++)
                    otf2[u, v] = OpticalTransferFunction.ModPow(otf[u, v]);

            for (int u = 0; u < height; u++)
                for (int v = 0; v < width; v++)
                    p[u, v] = p[u, v] * gamma;

            for (int u = 0; u < height; u++)
                for (int v = 0; v < width; v++)
                    otf2[u, v] = otf2[u, v] + p[u, v];

            for (int u = 0; u < height; u++)
                for (int v = 0; v < width; v++)
                    otf[u, v] = otfZ[u, v] / otf2[u, v];

            ConvolutionFilter cf = OpticalTransferFunction.Otf2psf(otf);

            return cf;
        }

 protected override IterateSinglePointResult? CreateSinglePoint(
     Complex current,
     Complex parameter,
     int maxIterations,
     double maxMagnitude)
 {
     return this.IterateSinglePoint(
         initial: current,
         offset: parameter,
         returnValue: current,
         maxIterations: maxIterations,
         maxMagnitude: maxMagnitude);
 }

        public void FourierBluesteinIsReversible(FourierOptions options)
        {
            var dft = new DiscreteFourierTransform();

            var samples = Generate.RandomComplex(0x7FFF, GetUniform(1));
            var work = new Complex[samples.Length];
            samples.CopyTo(work, 0);

            dft.BluesteinForward(work, options);
            Assert.IsFalse(work.ListAlmostEqual(samples, 6));

            dft.BluesteinInverse(work, options);
            AssertHelpers.ListAlmostEqual(samples, work, 10);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserCholesky"/> class. This object will compute the
        /// Cholesky factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not positive definite.</exception>
        public UserCholesky(Matrix<Complex> matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
            CholeskyFactor = matrix.Clone();
            var tmpColumn = new Complex[CholeskyFactor.RowCount];

            // Main loop - along the diagonal
            for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
            {
                // "Pivot" element
                var tmpVal = CholeskyFactor.At(ij, ij);

                if (tmpVal.Real > 0.0)
                {
                    tmpVal = tmpVal.SquareRoot();
                    CholeskyFactor.At(ij, ij, tmpVal);
                    tmpColumn[ij] = tmpVal;

                    // Calculate multipliers and copy to local column
                    // Current column, below the diagonal
                    for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                    {
                        CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
                        tmpColumn[i] = CholeskyFactor.At(i, ij);
                    }

                    // Remaining columns, below the diagonal
                    DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
                }
                else
                {
                    throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                }

                for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                {
                    CholeskyFactor.At(ij, i, Complex.Zero);
                }
            }
        }

        public static Vector<Complex> GenerateRandomDenseVector(int order)
        {
            // Fill a matrix with standard random numbers.
            var normal = new Distributions.Normal
                         {
                             RandomSource = new Random.MersenneTwister(1)
                         };
            var v = new DenseVector(order);
            for (var i = 0; i < order; i++)
            {
                v[i] = new Complex(normal.Sample(), normal.Sample());
            }

            // Generate a matrix which is positive definite.
            return v;
        }

        public void FourierDefaultTransformSatisfiesParsevalsTheorem(int count)
        {
            var samples = Generate.RandomComplex(count, GetUniform(1));

            var timeSpaceEnergy = (from s in samples select s.MagnitudeSquared()).Mean();

            var work = new Complex[samples.Length];
            samples.CopyTo(work, 0);

            // Default -> Symmetric Scaling
            Fourier.Forward(work);

            var frequencySpaceEnergy = (from s in work select s.MagnitudeSquared()).Mean();

            Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-12);
        }

        public void FourierDefaultTransformSatisfiesParsevalsTheorem(int count)
        {
            var samples = Sample.Random((u, v) => new Complex(u, v), _uniform, count);

            var timeSpaceEnergy = (from s in samples select s.MagnitudeSquared()).Mean();

            var work = new Complex[samples.Length];
            samples.CopyTo(work, 0);

            // Default -> Symmetric Scaling
            Transform.FourierForward(work);

            var frequencySpaceEnergy = (from s in work select s.MagnitudeSquared()).Mean();

            Assert.AreApproximatelyEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-12);
        }

 protected IterateSinglePointResult? IterateSinglePoint(
     Complex initial,
     Complex offset,
     Complex returnValue,
     int maxIterations,
     double maxMagnitude)
 {
     var z1 = initial;
     for (int i = 0; i < maxIterations; i++)
     {
         z1 = z1 * z1 + offset;
         if (z1.Magnitude > maxMagnitude)
         {
             return new IterateSinglePointResult { C = returnValue, Iterations = i };
         }
     }
     return null;
 }