8000 Add Determinant Gauss, Determinant - Kraut. · cp-algorithms/cp-algorithms@5501719 · GitHub
[go: up one dir, main page]
Skip to content

Commit 5501719

Browse files
author
Trung, Nguyen
committed
Add Determinant Gauss, Determinant - Kraut.
1 parent a9f533f commit 5501719

File tree

4 files changed

+156
-4
lines changed

4 files changed

+156
-4
lines changed

‎src/index.md‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,9 @@ especially popular in field of competitive programming.*
2424
- [Z-function](./string/z-function.html)
2525

2626
### Linear Algebra
27-
- [Gauss & System of linear equation](./algebra/linear-system-gauss.html)
27+
- [Gauss & System of linear equation](./linear_algebra/linear-system-gauss.html)
28+
- [Gauss & Determinant](./linear_algebra/determinant-gauss.html)
29+
- [Kraut & Determinant](./linear_algebra/determinant-kraut.html)
2830

2931
---
3032

‎src/linear_algebra/determinant-gauss.md‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,47 @@
1+
<!--?title Calculating Determinant of Matrix by Gauss -->
2+
3+
# Calculating the determinant of a matrix by Gauss
4+
5+
Problem: Given a matrix $A$ of size $N x N$. Compute its determinant.
6+
7+
## Algorithm
8+
9+
We use the ideas of [Gauss method for solving systems of linear equations](./linear_algebra/linear-system-gauss.html)
10+
11+
We will perform the same steps as in the solution of systems of linear equations, excluding only the division of the current line to $a_{ij}$. These operations will not change the absolute value of the determinant of the matrix. When we exchange two lines of the matrix, however, the sign of the determinant can change.
12+
13+
After applying Gauss on the matrix, we receive a diagonal matrix, whose determinant is just the product of the elements on the diagonal. The sign, as previously mentioned, can be determined by the number of exchanged rows (if odd, then the sign of the determinant should be reversed). Thus, we can use the Gauss algorithm to compute the determinant of hte matrix in comlexity $O(N^3)$.
14+
15+
It should be noted that if at some point, we do not find non-zero cell in current column, the algorithm should stop and returns 0.
16+
17+
## Implementation
18+
19+
```cpp
20+
const double EPS = 1E-9;
21+
int n;
22+
vector < vector<double> > a (n, vector<double> (n));
23+
24+
double det = 1;
25+
for (int i=0; i<n; ++i) {
26+
int k = i;
27+
for (int j=i+1; j<n; ++j)
28+
if (abs (a[j][i]) > abs (a[k][i]))
29+
k = j;
30+
if (abs (a[k][i]) < EPS) {
31+
det = 0;
32+
break;
33+
}
34+
swap (a[i], a[k]);
35+
if (i != k)
36+
det = -det;
37+
det *= a[i][i];
38+
for (int j=i+1; j<n; ++j)
39+
a[i][j] /= a[i][i];
40+
for (int j=0; j<n; ++j)
41+
if (j != i && abs (a[j][i]) > EPS)
42+
for (int k=i+1; k<n; ++k)
43+
a[j][k] -= a[i][k] * a[j][i];
44+
}
45+
46+
cout << det;
47+
```

‎src/linear_algebra/determinant-kraut.md‎

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
<!--?title Calculating the determinant using Kraut method-->
2+
3+
# Calculating the determinant using Kraut method in $O(N^3)$
4+
5+
In this article, we'll describe how to find the determinant of the matrix using Kraut method, which works in $O(N^3)$.
6+
7+
The Kraut algorithm finds decomposition of matrix $A$ as $A = L U$ where L is lower and and U is upper triangular matrix. Without loss of generality, we can assume that all the diagonal elements of L are equal to 1. But knowing these matrices, it is easy to calculate the determinant of A: it is equal to the product of all the elements on the main diagonal of the matrix U.
8+
9+
There is a theorem stating that any invertible matrix has a LU-decomposition, and it is unique, if and only if all its principle minors are non-zero. It should be recall that we consider only such decomposition in which the diagonal L consists of ones.
10+
11+
Let A be the matrix and N is its size. We will find the elements of the matrices L and U using the following steps:
12+
13+
1. Let $L_{i i} = 1$ for $i = 1, 2, ..., N$
14+
2. For each $j = 1, 2, ..., N$ perform:
15+
- For $i = 1, 2, ..., j$ find value $U_{ij}$:
16+
$U_{ij} = A_{ij} - Sum( L_{ik} * U_{kj} )$
17+
where the sum is over all $k = 1, 2, ..., i-1$.
18+
- Next, for $i = j+1, j+2, ..., N$ have:
19+
$L_{ij} = (A_{ij} - Sum( L_{ik} * U_{kj} )) / U_{jj}$,
20+
where the sum is taken over all $k = 1, 2, ..., j-1$.
21+
22+
## Implementation
23+
24+
```java
25+
static BigInteger det (BigDecimal a [][], int n) {
26+
try {
27+
28+
for (int i=0; i<n; i++) {
29+
boolean nonzero = false;
30+
for (int j=0; j<n; j++)
31+
if (a[i][j].compareTo (new BigDecimal (BigInteger.ZERO)) > 0)
32+
nonzero = true;
33+
if (!nonzero)
34+
return BigInteger.ZERO;
35+
}
36+
37+
BigDecimal scaling [] = new BigDecimal [n];
38+
for (int i=0; i<n; i++) {
39+
BigDecimal big = new BigDecimal (BigInteger.ZERO);
40+
for (int j=0; j<n; j++)
41+
if (a[i][j].abs().compareTo (big) > 0)
42+
big = a[i][j].abs();
43+
scaling[i] = (new BigDecimal (BigInteger.ONE)) .divide
44+
(big, 100, BigDecimal.ROUND_HALF_EVEN);
45+
}
46+
47+
int sign = 1;
48+
49+
for (int j=0; j<n; j++) {
50+
for (int i=0; i<j; i++) {
51+
BigDecimal sum = a[i][j];
52+
for (int k=0; k<i; k++)
53+
sum = sum.subtract (a[i][k].multiply (a[k][j]));
54+
a[i][j] = sum;
55+
}
56+
57+
BigDecimal big = new BigDecimal (BigInteger.ZERO);
58+
int imax = -1;
59+
for (int i=j; i<n; i++) {
60+
BigDecimal sum = a[i][j];
61+
for (int k=0; k<j; k++)
62+
sum = sum.subtract (a[i][k].multiply (a[k][j]));
63+
a[i][j] = sum;
64+
BigDecimal cur = sum.abs();
65+
cur = cur.multiply (scaling[i]);
66+
if (cur.compareTo (big) >= 0) {
67+
big = cur;
68+
imax = i;
69+
}
70+
}
71+
72+
if (j != imax) {
73+
for (int k=0; k<n; k++) {
74+
BigDecimal t = a[j][k];
75+
a[j][k] = a[imax][k];
76+
a[imax][k] = t;
77+
}
78+
79+
BigDecimal t = scaling[imax];
80+
scaling[imax] = scaling[j];
81+
scaling[j] = t;
82+
83+
sign = -sign;
84+
}
85+
86+
if (j != n-1)
87+
for (int i=j+1; i<n; i++)
88+
a[i][j] = a[i][j].divide
89+
(a[j][j], 100, BigDecimal.ROUND_HALF_EVEN);
90+
91+
}
92+
93+
BigDecimal result = new BigDecimal (1);
94+
if (sign == -1)
95+
result = result.negate();
96+
for (int i=0; i<n; i++)
97+
result = result.multiply (a[i][i]);
98+
99+
return result.divide
100+
(BigDecimal.valueOf(1), 0, BigDecimal.ROUND_HALF_EVEN).toBigInteger();
101+
}
102+
catch (Exception e) {
103+
return BigInteger.ZERO;
104+
}
105+
}
106+
```

src/algebra/linear-system-gauss.md renamed to src/linear_algebra/linear-system-gauss.md

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -189,6 +189,3 @@ Therefore, the resulting Gauss-Jordan solution must sometimes be improved by app
189189

190190
Thus, the solution turns into two-step: First, Gauss-Jordan algorithm is applied, and then a numerical method taking initial solution as solution in the first step.
191191

192-
## Translator's note
193-
194-
In my attempt of translating this article using Google Translate, there are some parts which was very hard to understand and translate properly. Thus in many places, I used my own words, and the translation is not exactly the same as original article. If you find any improvement, please contact me at ngthanhtrung23@gmail.com

0 commit comments

Comments
 (0)
0