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;        }    }}