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)