CP-templates
This documentation is automatically generated by online-judge-tools/verification-helper
icpc/eliminate_double.cpp
- View this file on GitHub
- Last update: 2025-03-22 15:11:28+08:00
Code
//source: KACTL
double det(vector<vector<double>>& a) {
int n = sz(a); double res = 1;
rep(i,0,n) {
int b = i;
rep(j,i+1,n) if (fabs(a[j][i]) > fabs(a[b][i])) b = j;
if (i != b) swap(a[i], a[b]), res *= -1;
res *= a[i][i];
if (res == 0) return 0;
rep(j,i+1,n) {
double v = a[j][i] / a[i][i];
if (v != 0) rep(k,i+1,n) a[j][k] -= v * a[i][k];
}
}
return res;
}
k
int matInv(vector<vector<double>>& A) {
int n = sz(A); vi col(n);
vector<vector<double>> tmp(n, vector<double>(n));
rep(i,0,n) tmp[i][i] = 1, col[i] = i;
rep(i,0,n) {
int r = i, c = i;
rep(j,i,n) rep(k,i,n)
if (fabs(A[j][k]) > fabs(A[r][c]))
r = j, c = k;
if (fabs(A[r][c]) < 1e-12) return i;
A[i].swap(A[r]); tmp[i].swap(tmp[r]);
rep(j,0,n)
swap(A[j][i], A[j][c]), swap(tmp[j][i], tmp[j][c]);
swap(col[i], col[c]);
double v = A[i][i];
rep(j,i+1,n) {
double f = A[j][i] / v;
A[j][i] = 0;
rep(k,i+1,n) A[j][k] -= f*A[i][k];
rep(k,0,n) tmp[j][k] -= f*tmp[i][k];
}
rep(j,i+1,n) A[i][j] /= v;
rep(j,0,n) tmp[i][j] /= v;
A[i][i] = 1;
}
/// forget A at this point, just eliminate tmp backward
for (int i = n-1; i > 0; --i) rep(j,0,i) {
double v = A[j][i];
rep(k,0,n) tmp[j][k] -= v*tmp[i][k];
}
rep(i,0,n) rep(j,0,n) A[col[i]][col[j]] = tmp[i][j];
return n; //rank
}
typedef vector<double> vd;
const double eps = 1e-12;
//Solves $A * x = b$. If there are multiple solutions, an arbitrary one is returned.
//Returns rank, or -1 if no solutions. Data in $A$ and $b$ is lost.
int solveLinear(vector<vd>& A, vd& b, vd& x) {
int n = sz(A), m = sz(x), rank = 0, br, bc;
if (n) assert(sz(A[0]) == m);
vi col(m); iota(all(col), 0);
rep(i,0,n) {
double v, bv = 0;
rep(r,i,n) rep(c,i,m)
if ((v = fabs(A[r][c])) > bv)
br = r, bc = c, bv = v;
if (bv <= eps) {
rep(j,i,n) if (fabs(b[j]) > eps) return -1;
break;
}
swap(A[i], A[br]);
swap(b[i], b[br]);
swap(col[i], col[bc]);
rep(j,0,n) swap(A[j][i], A[j][bc]);
bv = 1/A[i][i];
rep(j,i+1,n) {
double fac = A[j][i] * bv;
b[j] -= fac * b[i];
rep(k,i+1,m) A[j][k] -= fac*A[i][k];
}
rank++;
}
x.assign(m, 0);
for (int i = rank; i--;) {
b[i] /= A[i][i];
x[col[i]] = b[i];
rep(j,0,i) b[j] -= A[j][i] * b[i];
}
return rank; // (multiple solutions if rank < m)
}#line 1 "icpc/eliminate_double.cpp"
//source: KACTL
double det(vector<vector<double>>& a) {
int n = sz(a); double res = 1;
rep(i,0,n) {
int b = i;
rep(j,i+1,n) if (fabs(a[j][i]) > fabs(a[b][i])) b = j;
if (i != b) swap(a[i], a[b]), res *= -1;
res *= a[i][i];
if (res == 0) return 0;
rep(j,i+1,n) {
double v = a[j][i] / a[i][i];
if (v != 0) rep(k,i+1,n) a[j][k] -= v * a[i][k];
}
}
return res;
}
k
int matInv(vector<vector<double>>& A) {
int n = sz(A); vi col(n);
vector<vector<double>> tmp(n, vector<double>(n));
rep(i,0,n) tmp[i][i] = 1, col[i] = i;
rep(i,0,n) {
int r = i, c = i;
rep(j,i,n) rep(k,i,n)
if (fabs(A[j][k]) > fabs(A[r][c]))
r = j, c = k;
if (fabs(A[r][c]) < 1e-12) return i;
A[i].swap(A[r]); tmp[i].swap(tmp[r]);
rep(j,0,n)
swap(A[j][i], A[j][c]), swap(tmp[j][i], tmp[j][c]);
swap(col[i], col[c]);
double v = A[i][i];
rep(j,i+1,n) {
double f = A[j][i] / v;
A[j][i] = 0;
rep(k,i+1,n) A[j][k] -= f*A[i][k];
rep(k,0,n) tmp[j][k] -= f*tmp[i][k];
}
rep(j,i+1,n) A[i][j] /= v;
rep(j,0,n) tmp[i][j] /= v;
A[i][i] = 1;
}
/// forget A at this point, just eliminate tmp backward
for (int i = n-1; i > 0; --i) rep(j,0,i) {
double v = A[j][i];
rep(k,0,n) tmp[j][k] -= v*tmp[i][k];
}
rep(i,0,n) rep(j,0,n) A[col[i]][col[j]] = tmp[i][j];
return n; //rank
}
typedef vector<double> vd;
const double eps = 1e-12;
//Solves $A * x = b$. If there are multiple solutions, an arbitrary one is returned.
//Returns rank, or -1 if no solutions. Data in $A$ and $b$ is lost.
int solveLinear(vector<vd>& A, vd& b, vd& x) {
int n = sz(A), m = sz(x), rank = 0, br, bc;
if (n) assert(sz(A[0]) == m);
vi col(m); iota(all(col), 0);
rep(i,0,n) {
double v, bv = 0;
rep(r,i,n) rep(c,i,m)
if ((v = fabs(A[r][c])) > bv)
br = r, bc = c, bv = v;
if (bv <= eps) {
rep(j,i,n) if (fabs(b[j]) > eps) return -1;
break;
}
swap(A[i], A[br]);
swap(b[i], b[br]);
swap(col[i], col[bc]);
rep(j,0,n) swap(A[j][i], A[j][bc]);
bv = 1/A[i][i];
rep(j,i+1,n) {
double fac = A[j][i] * bv;
b[j] -= fac * b[i];
rep(k,i+1,m) A[j][k] -= fac*A[i][k];
}
rank++;
}
x.assign(m, 0);
for (int i = rank; i--;) {
b[i] /= A[i][i];
x[col[i]] = b[i];
rep(j,0,i) b[j] -= A[j][i] * b[i];
}
return rank; // (multiple solutions if rank < m)
}