using System;

namespace Zhou.CSharp.Algorithm
{
    /// <summary>
    /// 矩阵类
    /// 作者：周长发
    /// 改进：深度混淆
    /// https://blog.csdn.net/beijinghorn
    /// </summary>
    public partial class Matrix
    {

        /// <summary>
        /// 求实对称矩阵特征值与特征向量的雅可比过关法
        /// </summary>
        /// <param name="src">源矩阵</param>
        /// <param name="dblEigenValue">一维数组，长度为矩阵的阶数，返回时存放特征值</param>
        /// <param name="mtxEigenVector">返回时存放特征向量矩阵，其中第i列为与数组dblEigenValue中第j个特征值对应的特征向量</param>
        /// <param name="eps">计算精度</param>
        /// <returns>求解是否成功</returns>
        public static bool ComputeEvJacobi(Matrix src, out double[] dblEigenValue, out Matrix mtxEigenVector, double eps = 1.0E-7)
        {
            int p, q, u, w, t, s;
            double ff, fm, cn, sn, omega, x, y, d;

            int n = src.Columns;
            dblEigenValue = new double[n];
            mtxEigenVector = new Matrix(n, n);

            for (int i = 0; i < n; i++)
            {
                mtxEigenVector[i * n + i] = 1.0;
                for (int j = 0; j < n; j++)
                {
                    if (i != j)
                    {
                        mtxEigenVector[i * n + j] = 0.0;
                    }
                }
            }
            ff = 0.0;
            for (int i = 1; i < n; i++)
            {
                for (int j = 0; j < i; j++)
                {
                    d = src[i * n + j];
                    ff = ff + d * d;
                }
            }
            if (Math.Abs(2.0 * ff) < float.Epsilon)
            {
                return false;
            }
            ff = Math.Sqrt(2.0 * ff);
            ff = ff / (1.0 * n);

            bool nextLoop = false;
            while (true)
            {
                for (int i = 1; i < n; i++)
                {
                    for (int j = 0; j < i; j++)
                    {
                        d = Math.Abs(src[i * n + j]);
                        if (d > ff)
                        {
                            p = i;
                            q = j;

                            u = p * n + q;
                            w = p * n + p;
                            t = q * n + p;
                            s = q * n + q;
                            x = -src[u];
                            y = (src[s] - src[w]) / 2.0;
                            if (Math.Abs(x) < float.Epsilon || Math.Abs(y) < float.Epsilon)
                            {
                                return false;
                            }
                            omega = x / Math.Sqrt(x * x + y * y);
                            if (y < 0.0)
                            {
                                omega = -omega;
                            }
                            if (Math.Abs(omega - 1.0) < float.Epsilon)
                            {
                                return false;
                            }
                            sn = 1.0 + Math.Sqrt(1.0 - omega * omega);
                            sn = omega / Math.Sqrt(2.0 * sn);
                            if (Math.Abs(sn - 1.0) < float.Epsilon)
                            {
                                return false;
                            }
                            cn = Math.Sqrt(1.0 - sn * sn);
                            fm = src[w];
                            src[w] = fm * cn * cn + src[s] * sn * sn + src[u] * omega;
                            src[s] = fm * sn * sn + src[s] * cn * cn - src[u] * omega;
                            src[u] = 0.0;
                            src[t] = 0.0;

                            for (int jj = 0; jj < n; jj++)
                            {
                                if ((jj != p) && (jj != q))
                                {
                                    u = p * n + jj;
                                    w = q * n + jj;
                                    fm = src[u];
                                    src[u] = fm * cn + src[w] * sn;
                                    src[w] = -fm * sn + src[w] * cn;
                                }
                            }

                            for (int ii = 0; ii < n; ii++)
                            {
                                if ((ii != p) && (ii != q))
                                {
                                    u = ii * n + p;
                                    w = ii * n + q;
                                    fm = src[u];
                                    src[u] = fm * cn + src[w] * sn;
                                    src[w] = -fm * sn + src[w] * cn;
                                }
                            }

                            for (int ii = 0; ii < n; ii++)
                            {
                                u = ii * n + p;
                                w = ii * n + q;
                                fm = mtxEigenVector[u];
                                mtxEigenVector[u] = fm * cn + mtxEigenVector[w] * sn;
                                mtxEigenVector[w] = -fm * sn + mtxEigenVector[w] * cn;
                            }

                            nextLoop = true;
                            break;
                        }
                    }

                    if (nextLoop)
                    {
                        break;
                    }
                }
                if (nextLoop)
                {
                    nextLoop = false;
                    continue;
                }

                nextLoop = false;

                // 如果达到精度要求，退出循环，返回结果
                if (ff < eps)
                {
                    for (int i = 0; i < n; ++i)
                    {
                        dblEigenValue[i] = src[i, i];
                    }
                    return true;
                }

                ff = ff / (1.0 * n);
            }
        }
    }
}