MyMathLib系列(行列式计算2)

/// <summary>
    /// 行列式计算,本程序属于MyMathLib的一部分。欢迎使用，參考，提意见。
    /// 有时间用函数语言改写，做自己得MathLib,里面的算法经过验证，但没经过
    /// 严格測试，如需參考，请谨慎.
    /// </summary>
    public static partial class LinearAlgebra
    { /// <summary>
        /// 获取指定i,j的余子式
        /// </summary>
        /// <param name="Determinants">N阶行列式</param>
        /// <param name="i">第i行</param>
        /// <param name="j">第j列</param>
        /// <returns>计算结果</returns>
        public static T[,] GetDeterminantMij<T>(T[,] Determinants, int i, int j)
        {
            var theN = Determinants.GetLength(0);
            var theNewDeter = new T[theN - 1, theN - 1];
            int theI = -1;

            for (int k = 0; k < theN; k++)
            {
                if (k == i - 1)
                {
                    continue;
                }
                theI++;
                int theJ = -1;
                for (int l = 0; l < theN; l++)
                {
                    if (l == j - 1)
                    {
                        continue;
                    }
                    theJ++;
                    theNewDeter[theI, theJ] = Determinants[k, l];
                }
            }
            return theNewDeter;
        }
        /// <summary>
        /// 获取指定i,j的余子式
        /// </summary>
        /// <param name="Determinants">N阶行列式</param>
        /// <param name="Rows">要取得行</param>
        /// <param name="Cols">要取得列</param>
        /// <returns>计算结果</returns>
        public static T[,] GetDeterminantMij<T>(T[,] Determinants, int[] Rows, int[] Cols)
        {
            if (Rows.Length != Cols.Length)
            {
                throw new Exception("所取行数和列数必须相等!");
            }
            var theN = Determinants.GetLength(0);
            var theNewN = theN - Rows.Length;
            var theNewDeter = new T[theNewN, theNewN];
            int theI = -1;

            for (int k = 0; k < theN; k++)
            {
                if (Rows.Contains(k + 1))
                {
                    continue;
                }
                theI++;
                int theJ = -1;
                for (int l = 0; l < theN; l++)
                {
                    if (Cols.Contains(l + 1))
                    {
                        continue;
                    }
                    theJ++;
                    theNewDeter[theI, theJ] = Determinants[k, l];
                }
            }
            return theNewDeter;
        }
        /// <summary>
        /// 获取指定k阶子式N
        /// </summary>
        /// <param name="Determinants">N阶行列式</param>
        /// <param name="Rows">要取得行</param>
        /// <param name="Cols">要取得列</param>
        /// <returns>计算结果</returns>
        public static T[,] GetDeterminantKN<T>(T[,] Determinants, int[] Rows, int[] Cols)
        {
            if (Rows.Length != Cols.Length)
            {
                throw new Exception("所取行数和列数必须相等!");
            }
            var theNewN = Rows.Length;
            var theNewDeter = new T[theNewN, theNewN];
            for (int k = 0; k < Rows.Length; k++)
            {
                for (int l = 0; l < Cols.Length; l++)
                {
                    theNewDeter[k, l] = Determinants[Rows[k] - 1, Cols[l] - 1];
                }
            }
            return theNewDeter;
        }
        /// <summary>
        /// 计算余子式的符号。

        /// </summary>
        /// <param name="i"></param>
        /// <param name="j"></param>
        /// <returns></returns>
        public static int CalcDeterMijSign(int i, int j)
        {
            int theSign = 1;
            if ((i + j) % 2 == 1)
            {
                theSign = -1;
            }
            return theSign;
        }
        /// <summary>
        /// 计算余子式的符号。
        /// </summary>
        /// <param name="i"></param>
        /// <param name="j"></param>
        /// <returns></returns>
        public static int CalcDeterMijSign(int[] Rows, int[] Cols)
        {
            int theSign = 1;
            var theSum = Rows.Sum() + Cols.Sum();
            if (theSum % 2 == 1)
            {
                theSign = -1;
            }
            return theSign;
        }
        /// <summary>
        /// 降阶法计算行列式
        /// </summary>
        /// <param name="Determinants">N阶行列式</param>
        /// <param name="ZeroOptimization">是否0优化</param>
        /// <returns>计算结果</returns>
        public static decimal CalcDeterminantAij(decimal[,] Determinants, bool ZeroOptimization = false)
        {
            var theN = Determinants.GetLength(0);
            //假设为2阶，直接计算
            if (theN == 2)
            {
                return Determinants[0, 0] * Determinants[1, 1] - Determinants[0, 1] * Determinants[1, 0];
            }
            if (ZeroOptimization)
            {
                //找0最多的行
                int theRowIndex = 0;
                int theMaxZeroCountR = -1;
                for (int i = 0; i < theN; i++)
                {
                    int theZeroNum = 0;
                    for (int j = 0; j < theN; j++)
                    {
                        if (Determinants[i, j] == 0)
                        {
                            theZeroNum++;
                        }
                    }
                    if (theZeroNum > theMaxZeroCountR)
                    {
                        theRowIndex = i;
                        theMaxZeroCountR = theZeroNum;
                    }
                }
                //找0最多的列
                int theColIndex = 0;
                int theMaxZeroCountC = -1;
                for (int i = 0; i < theN; i++)
                {
                    int theZeroNum = 0;
                    for (int j = 0; j < theN; j++)
                    {
                        if (Determinants[j, i] == 0)
                        {
                            theZeroNum++;
                        }
                    }
                    if (theZeroNum > theMaxZeroCountC)
                    {
                        theColIndex = i;
                        theMaxZeroCountC = theZeroNum;
                    }
                }
                if (theMaxZeroCountR >= theMaxZeroCountC)
                {
                    decimal theRetDec = 0;
                    //第i=theRowIndex+1行展开
                    int i = theRowIndex + 1;
                    for (int j = 1; j <= theN; j++)
                    {
                        var theSign = CalcDeterMijSign(i, j);
                        var theNewMij = GetDeterminantMij(Determinants, i, j);
                        theRetDec += theSign * Determinants[i - 1, j - 1] * CalcDeterminantAij(theNewMij, ZeroOptimization);
                    }
                    return theRetDec;
                }
                else
                {
                    decimal theRetDec = 0;
                    //第j=theColIndex+1列展开
                    int j = theColIndex + 1;
                    for (int i = 1; i <= theN; i++)
                    {
                        var theSign = CalcDeterMijSign(i, j);
                        var theNewMij = GetDeterminantMij(Determinants, i, j);
                        theRetDec += theSign * Determinants[i, j] * CalcDeterminantAij(theNewMij, ZeroOptimization);
                    }
                    return theRetDec;
                }
            }
            else
            {
                //採用随机法展开一行
                var i = new Random().Next(1, theN);
                decimal theRetDec = 0;
                for (int j = 1; j <= theN; j++)
                {
                    var theSign = CalcDeterMijSign(i, j);
                    var theNewMij = GetDeterminantMij(Determinants, i, j);
                    theRetDec += theSign * Determinants[i-1, j-1] * CalcDeterminantAij(theNewMij, ZeroOptimization);
                }
                return theRetDec;
            }
        }

        /// <summary>
        /// 计算范德蒙行列式
        /// </summary>
        /// <param name="Determinants">范德蒙行列式简记序列</param>
        /// <returns>计算结果</returns>
        public static decimal CalcVanDerModeDeter(decimal[] VanDerModeDeter)
        {
            var theN = VanDerModeDeter.Length;
            if (theN == 1)
            {
                return 1;
            }
            decimal theRetDec = 1;
            for (int i = 0; i < theN; i++)
            {
                for (int j = i + 1; j < theN; j++)
                {
                    theRetDec *= (VanDerModeDeter[j] - VanDerModeDeter[i]);
                }
            }
            return theRetDec;
        }
        /// <summary>
        /// 获取奇数序列
        /// </summary>
        /// <param name="N"></param>
        /// <returns></returns>
        private static int[] GetLaplaceRowsOdd(int N)
        {
            var theRet = new List<int>();
            for (int i = 0; i < N; i = i + 2)
            {
                theRet.Add(i + 1);
            }
            return theRet.ToArray();
        }

        /// <summary>
        /// 依据拉普拉斯定理计算行列式值。

        /// </summary>
        /// <param name="Determinants">N阶行列式</param>
        /// <param name="Rows">初始展开行,里面採用奇数行展开</param>
        /// <returns>计算结果</returns>
        public static decimal CalcDeterByLaplaceLaw(decimal[,] Determinants, int[] Rows)
        {
            var n = Determinants.GetLength(0);
            var k = Rows.Length;
            //假设阶数小于3，则不是必需採用拉普拉斯展开
            if (n <= 3)
            {
                return CalcDeterminantAij(Determinants, false);
            }
            //从P（theN,theK)
            var theRetList = GetCombination(n, k);
            decimal theRetDec = 0;
            foreach (var theCols in theRetList)
            {
                var theSign = CalcDeterMijSign(Rows, theCols.ToArray());
                var theKN = GetDeterminantKN(Determinants, Rows, theCols.ToArray());
                var theN = GetDeterminantMij(Determinants, Rows, theCols.ToArray());
                decimal theRetKN = 0;
                //假设剩余阶数>4则採用随机半数处理.
                if (n - k >= 4)
                {
                    var theRows = GetLaplaceRowsOdd(n - k);
                    theRetKN = CalcDeterByLaplaceLaw(theKN, theRows);
                }
                else
                {
                    theRetKN = CalcDeterminantAij(theKN);
                }
                decimal theRetAk = 0;
                if (k >= 4)
                {
                    var theRows = GetLaplaceRowsOdd(k);
                    theRetAk = CalcDeterByLaplaceLaw(theN, theRows);
                }
                else
                {
                    theRetAk = CalcDeterminantAij(theN);
                }
                theRetDec += theSign * theRetKN * theRetAk;
            }
            return theRetDec;
        }
        /// <summary>
        /// 从N个数中取k个数的组合结果。考虑到组合数没有顺序区分，因此仅仅要考虑从小
        /// 到大的排列下的组合情况就可以，另外，假设组合也不用考虑元素反复的
        /// 问题。假设有反复数，仅仅要除重就可以。

        /// </summary>
        /// <param name="N">N个数1-N</param>
        /// <param name="k">取K个</param>
        /// <returns></returns>
        public static List<List<int>> GetCombination(int N, int k)
        {
            var theList = new List<int>();
            for (int i = 1; i <= N; i++)
            {
                theList.Add(i);
            }
            return GetCombination(theList, k);
        }
        /// <summary>
        /// 从N个中取k个数,算法原理C(N,k)=C(N-1,k)+ (a + C(Na-1,k-1));当中Na是N中去掉a后的集合.
        /// </summary>
        /// <param name="N">元素总个数</param>
        /// <param name="k">取k个</param>
        /// <returns></returns>
        public static List<List<int>> GetCombination(List<int> N, int k)
        {
            if (k==0)
            {
                return null;
            }
            if (N.Count < k)
            {
                return null;
            }
            if (k == 1)
            {
                var theResultsList = new List<List<int>>();
                foreach (var theN in N)
                {
                    var theList = new List<int>();
                    theList.Add(theN);
                    theResultsList.Add(theList);
                }
                return theResultsList;
            }
            if (N.Count == k)
            {
                var theResultsList = new List<List<int>>();
                var theList = new List<int>();
                theList.AddRange(N);
                theResultsList.Add(theList);
                return theResultsList;
            }
            var theRet3 = new List<List<int>>();
            int theLeft = N[0];
            var theRight = new List<int>();
            theRight.AddRange(N);
            theRight.Remove(N[0]);

            var theRet2 = GetCombination(theRight, k);
            theRet3.AddRange(theRet2);

            theRet2 = GetCombination(theRight, k - 1);
            for (int n = 0; n < theRet2.Count; n++)
            {
                var theList = new List<int>();
                theList.Add(theLeft);
                theList.AddRange(theRet2[n]);
                theRet3.Add(theList);
            }
            return theRet3;
        }
    }
}