行列式の計算方法は 掃き出し法で連立一次方程式を解く 逆行列と行列式を求めるで解説しています。ここではある行の定数倍を別の行に加えても行列式は同じであるという性質と余因子展開を用いて行数を減らすことで行列式を求めています。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 |
double Determinant(double[,] matrix) { if (matrix.GetLength(0) != matrix.GetLength(1)) throw new Exception("正方行列ではない！！"); int n = matrix.GetLength(0); double[,] temp = new double[n, n]; for (int row = 0; row < n; row++) { for (int col = 0; col < n; col++) temp[row, col] = matrix[row, col]; } double scale = 1; while (true) { int rowCount = temp.GetLength(0); int colCount = temp.GetLength(1); if (rowCount == 1) return temp[0, 0] * scale; double leftTop = temp[0, 0]; if (leftTop != 0) scale *= leftTop; else { for (int row = 0; row < rowCount; row++) { if (temp[row, 0] != 0) { for (int col = 0; col < colCount; col++) { double d = temp[0, col]; temp[0, col] = temp[row, col]; temp[row, col] = d; } break; } } leftTop = temp[0, 0]; if (leftTop != 0) scale *= -leftTop; else return 0; } for (int col = 0; col < colCount; col++) temp[0, col] /= leftTop; for (int row = 1; row < rowCount; row++) { double d = temp[row, 0]; for (int col = 0; col < colCount; col++) temp[row, col] -= temp[0, col] * d; } double[,] submatrix = new double[rowCount - 1, colCount - 1]; for (int row = 1; row < rowCount; row++) { for (int col = 1; col < colCount; col++) submatrix[row - 1, col - 1] = temp[row, col]; } temp = submatrix; } } |
すべての要素が整数であるなら行列式も整数になりますが、この処理では割り算を使うため、誤差が出る可能性があります。また要素が整数である行列の行列式を 10^9 + 7で割ったときの余りを求めるときはどうすればよいのでしょうか？
逆元
A[0] ～ A[n] の合計を m で割ったときの剰余を考える時、(A[0]＋A[1]＋A[2]) % m = ((A[0]＋A[1]) % m＋A[2]) % m なので A[0] ～ A[n] の合計を計算しなくてもループのなかで以下のような処理をすれば解を得ることができます。掛け算も同様です。これによってオーバーフローを気にする必要はなくなります。
|
1 2 3 4 5 6 |
long ans = 0; foreach (int a in A) { ans += a; ans %= m; } |
引き算のときもほぼ同様ですが注意が必要です。
1 % 4 は 1、-3 は 1 から 4 を引いたものなので 1 % 4 == -3 % 4 となってほしいのですが、そうはなりません。-3 % 4 == -3 です。なので演算子 % を使用したときは結果が負数のときは m を足さなければなりません。
|
1 2 3 4 5 6 7 8 |
long ans = 0; foreach (int a in A) { ans -= a; ans %= m; if(ans < 0) ans += m; } |
問題は割り算です。割り算のときはそうはいきません。
ところで割り算は逆数を掛けることと同じです。剰余を考えるときも逆数のようなものが存在します。
a×bを m で割ったときの余りが 1 になるとき b を a の逆元といいます。
例えば mod 13 において 2 の逆元は 7 です。（2×7を13で割ったときの余りは1なので）
逆元を求めるには？
逆元を求めるには拡張ユークリッドの互助法を用います。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
long Inverse(long a, long mod) { if (a < 0) a += mod; long b = mod, u = 1, v = 0; while (b > 0) { long t = a / b; a -= t * b; swap(ref a, ref b); u -= t * v; swap(ref u, ref v); } u %= mod; if (u < 0) u += mod; return u; void swap(ref long a, ref long b) { long c = a; a = b; b = c; } } |
行列式を計算する
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 |
int Determinant(int[,] mat, int mod) { if (mat.GetLength(0) != mat.GetLength(1)) throw new Exception("正方行列ではない！！"); void ModRef(ref long val) { val %= mod; if (val < 0) val += mod; } int n = mat.GetLength(0); long[,] temp = new long[n, n]; for (int row = 0; row < n; row++) { for (int col = 0; col < n; col++) { temp[row, col] = mat[row, col]; ModRef(ref temp[row, col]); } } long scale = 1; while (true) { int rowCount = temp.GetLength(0); int colCount = temp.GetLength(1); if (rowCount == 1) { temp[0, 0] *= scale; ModRef(ref temp[0, 0]); return (int)temp[0, 0]; } ModRef(ref temp[0, 0]); long leftTop = temp[0, 0]; if (leftTop != 0) { scale *= leftTop; ModRef(ref scale); } else { for (int row = 0; row < rowCount; row++) { if (temp[row, 0] != 0) { for (int col = 0; col < colCount; col++) { long t = temp[0, col]; temp[0, col] = temp[row, col]; temp[row, col] = t; } break; } } leftTop = temp[0, 0]; if (leftTop != 0) { long t = -leftTop; ModRef(ref t); scale *= t; ModRef(ref scale); } else return 0; } for (int col = 0; col < colCount; col++) { temp[0, col] *= Inverse(leftTop, mod); ModRef(ref temp[0, col]); } for (int row = 1; row < rowCount; row++) { long t1 = temp[row, 0]; for (int col = 0; col < colCount; col++) { long t2 = -temp[0, col]; ModRef(ref t2); temp[row, col] += t1 * t2; ModRef(ref temp[row, col]); } } long[,] submatrix = new long[rowCount - 1, colCount - 1]; for (int row = 1; row < rowCount; row++) { for (int col = 1; col < colCount; col++) submatrix[row - 1, col - 1] = temp[row, col]; } temp = submatrix; } } |
使ってみる
問題の趣旨は以下のとおりです。
2次元配列 A が与えられる。
N 頂点のグラフ があり、A[i, j] == 1 のとき i と頂点 j を直接結ぶ辺が存在する。
A[i, j] == 0 のときは存在してはならない。
A[i, j] == -1 のときはどちらでもよい。
この条件で全域木は何通り作れるだろうか？
ラプラシアン行列と行列木定理
各行、列がそれぞれ頂点に対応した行列で、対角成分にはその頂点の次数、非対角成分については枝がある部分に ?1、ない部分に 0 を格納した正方行列をラプラシアン行列と言います。
そして行列木定理では「無向グラフ G の全域木の個数はグラフのラプラシアン行列の任意の余因子と等しい」とされています。
ラプラシアン行列の生成
そこで与えられた2次元配列からラプラシアン行列を生成します。そのまえに必ず辺が存在する部分だけで閉路が存在する場合、絶対に木にはならないのでその確認をおこないます。
そのあとそれぞれの連結成分をひとつの頂点に圧縮して、頂点と頂点の間に辺を張って全域木をつくる問題にするのですが、まず2次元配列において辺が存在してもしなくてもかまわない頂点のうち、連結成分の代表元同士で張ることにします。そのあと連結成分の内部で代表元とそれ以外の頂点で辺を張ります。こうすることで題意を満たしたラプラシアン行列を得ることができます。
連結成分を調べるにあたって ac-library-csharp を使っています。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 |
class Program { static void Main() { int N = int.Parse(Console.ReadLine()); int[,] A = new int[N, N]; for (int r = 0; r < N; r++) { int[] vs = Console.ReadLine().Split().Select(_ => long.Parse(_)).ToArray(); for (int c = 0; c < N; c++) A[r, c] = vs[c]; } AtCoder.Dsu dsu = new AtCoder.Dsu(N); for (int a = 0; a < N; a++) { for (int b = a + 1; b < N; b++) { if (A[a, b] == 1) { if(dsu.Same(a, b)) // 閉路があるので木はできない { Console.WriteLine(0); return; } dsu.Merge(a, b); } } } int[,] m = new int[N, N]; // 辺を張るが代表元同士でおこなう for (int a = 0; a < N; a++) { for (int b = a + 1; b < N; b++) { if (A[a, b] == -1) { int al = dsu.Leader(a); int bl = dsu.Leader(b); m[al, al]++; // ラプラシアン行列の対角成分にはその頂点の次数 m[bl, bl]++; m[al, bl]--; // それ以外の成分は枝がある部分に -1 m[bl, al]--; } } } // 連結成分内で代表元とその他の元で辺を張りあう for (int i = 0; i < N; i++) { int leader = dsu.Leader(i); m[leader, leader]++; m[i, i]++; m[leader, i]--; m[i, leader]--; } int[,] m2 = new int[N - 1, N - 1]; for (int row = 0; row < N - 1; row++) { for (int col = 0; col < N - 1; col++) m2[row, col] = m[row, col]; } int mod = (int)Math.Pow(10, 9) + 7; Console.WriteLine(Determinant(m2, mod)); } } |
