首页 > 代码库 > MyMathLib系列(行列式计算2)
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;
}
}
}
MyMathLib系列(行列式计算2)