一套矩阵库，主要实现功能有求特征值&特征向量，求合同标准型。
定义了实矩阵，多项式，多项式矩阵，向量四个类
注释比较详尽
/* Copyright © 2016 by dhd All Right Reserved. */ #include<bits/stdc++.h> #define rep(i, j, k) for(int i = j;i <= k;i++) #define repm(i, j, k) for(int i = j;i >= k;i--) #define mem(a) memset(a, 0, sizeof(a)) using namespace std; const double eps = 1e-8; const double pi = acos(-1); const double EPS = 1e-12; const double inf = 1e12; inline int sign(double x) { return x < -EPS ? -1 : x > EPS; } //get the value of a polynomial(the variable equals to x) inline double get(const vector<double> &coef, double x) { double e = 1 , s = 0; rep(i, 0, coef.size() - 1) s += e * coef[i], e *= x; return s; } //get the equation's root between lo & hi double find(const vector<double> &coef, int n, double lo, double hi) { int sign_lo, sign_hi; if((sign_lo = sign(get(coef, lo))) == 0) return lo; if((sign_hi = sign(get(coef, hi))) == 0) return hi; if(sign_lo * sign_hi > 0) return inf; for(int step = 0; step < 100 && hi - lo > EPS; ++step) { double mid = (lo + hi) * 0.5; int sign_mid = sign(get(coef, mid)); if(sign_mid == 0) return mid; if(sign_mid * sign_lo > 0) lo = mid; else hi = mid; } return (lo + hi) * 0.5; } //get all roots of the equation vector<double> solve(vector<double> coef, int n) { vector<double> ret; if(n == 1) { if(sign(coef[1])) ret.push_back(-coef[0] / coef[1]); return ret; } vector<double> dcoef(n); rep(i, 0, n-1) dcoef[i] = coef[i+1] * (i + 1); vector<double> droot = solve(dcoef, n - 1); droot.insert(droot.begin(), -inf); droot.push_back(inf); for(int i = 0; i + 1 < droot.size(); ++i) { double tmp = find(coef, n, droot[i], droot[i + 1]); if(tmp < inf) ret.push_back(tmp); } return ret; } //a class of polynomial struct poly{ int n; double num[210]; void clear() { rep(i, 0, 209) num[i] = 0; } poly() {} //using a vactor to generate a polynomial poly(vector<double> vec): n(vec.size() - 1) { clear(); rep(i, 0, vec.size() - 1) num[i] = vec[i]; } void print() { cout << endl; cout << num[n] << "x^" << n; repm(i, n-1, 0) { if(num[i] >= 0) cout << "+" << num[i] << "x^" << i; else cout << num[i] << "x^" << i; } cout << endl; } //return the addition of two polynomials friend poly operator + (poly a, poly b) { poly ans; ans.n = max(a.n, b.n); ans.clear(); rep(i, 0, ans.n) ans.num[i] = a.num[i] + b.num[i]; while(abs(ans.num[ans.n]) < eps && ans.n > 0) ans.n--; return ans; } //return the product of a polynomial and a number friend poly operator * (poly a, double b) { poly ans = a; rep(i, 0, ans.n) ans.num[i] *= b; ans.n = 200; while(abs(ans.num[ans.n]) < eps && ans.n > 0) ans.n--; return ans; } //return the subtraction of two polynomials friend poly operator - (poly a, poly b) { return a + (b * -1.0); } //return the product of two polynomials friend poly operator * (poly a, poly b) { poly ans; ans.n = 200; ans.clear(); rep(i, 0, a.n) rep(j, 0, b.n) ans.num[i + j] += a.num[i] * b.num[j]; while(abs(ans.num[ans.n]) < eps && ans.n > 0) ans.n--; return ans; } //return the ratio of two polynomials friend poly operator / (poly a, poly b) { poly ans; ans.n = 200; ans.clear(); while(a.n || abs(a.num[0]) > eps) { int dif = a.n - b.n; vector<double> vec; rep(i, 1, dif) vec.push_back(0); vec.push_back(a.num[a.n] / b.num[b.n]); poly part = poly(vec); a = a - (part * b); ans = ans + part; } while(abs(ans.num[ans.n]) < eps&&ans.n > 0) ans.n--; return ans; } bool empty() { if(n == 0 && abs(num[0]) < eps) return 1; return 0; } //return if a can be devided by b friend bool can_be_devided (poly a, poly b) { while(a.n >= b.n && !a.empty()) { int dif = a.n - b.n; vector<double> vec; rep(i, 1, dif) vec.push_back(0); vec.push_back(a.num[a.n] / b.num[b.n]); poly part = poly(vec); a = a - (part * b); } if(!a.empty()) return 0; return 1; } //get the roots of the equation(polynomial==0) vector<double> solu() { vector<double> vec; rep(i, 0, n) vec.push_back(num[i]); return solve(vec, n); } }; //a class of matrix; the elements are reals struct matrix{ int n, m; double num[101][101]; //scan the arguments & elements void scan(){ cin >> n >> m; rep(i, 1, n) rep(j, 1, m) cin >> num[i][j]; } //fill the matrix with element 0 void clear(){ rep(i, 1, 100) rep(j, 1, 100) num[i][j] = 0; } //let the matrix equals to I(or E) void normalize(){ clear(); rep(i, 1, n) num[i][i] = 1.0; } matrix () { clear(); } //print the arguments and elements void print(){ cout << endl; if(n == 0 || m == 0) { cout << "matrix is empty" << endl; return; } cout << n << " " << m << endl; rep(i, 1, n) { rep(j, 1, m) if(abs(num[i][j]) > eps) cout << num[i][j] << " "; else cout << "0 "; cout << endl; } } //swap the two lines of the matrix void swap_hang(int i, int j) { rep(cnt, 1, m) swap(num[i][cnt], num[j][cnt]); } //swap the two rows of the matrix void swap_lie(int i, int j) { rep(cnt, 1, n) swap(num[cnt][i], num[cnt][j]); } //update the j(th) line with the value of line(j)-mul*line(i) void minus_hang(int i, int j, double mul) { rep(cnt, 1, m) num[j][cnt] -= num[i][cnt] * mul; } //update the j(th) row with the value of row(j)-mul*row(i) void minus_lie(int i, int j, double mul){ rep(cnt, 1, n) num[cnt][j] -= num[cnt][i] * mul; } //update the i(th) line with the value of mul*line(i) void mul_hang(int i, double mul) { rep(cnt, 1, m) num[i][cnt] *= mul; } //update the i(th) row with the value of mul*row(i) void mul_lie(int i, double mul) { rep(cnt, 1, n) num[cnt][i] *= mul; } //return the addition of two matrixs matrix operator + (matrix b){ matrix x = (*this); if(n == b.n && m == b.m) rep(i, 1, n) rep(j, 1, m) x.num[i][j] += b.num[i][j]; cout << "can't add, will return the first matrix" << endl; return x; } //return the subtraction of two matrixs matrix operator - (matrix b){ matrix x = (*this); if(n == b.n && m == b.m) rep(i, 1, n) rep(j, 1, m) x.num[i][j] -= b.num[i][j]; else cout << "can't subtract, will return the first matrix" << endl; return x; } //return the product of two matrixs matrix operator * (matrix b){ matrix x = (*this); if(m == b.n) { x.n = n; x.m = b.m; x.clear(); rep(i, 1, x.n) rep(j, 1, x.m) rep(k, 1, m) x.num[i][k] += num[i][j] * b.num[j][k]; } else cout << "can't multiply, will return the first matrix" << endl; return x; } //return the product of a matrix and a number friend matrix operator *(matrix a, double b) { rep(i, 1, a.n) rep(j, 1, a.m) a.num[i][j] *= b; return a; } //return the transposition matrix of the original matrix matrix transposition(){ matrix a; a.m = n, a.n = m; rep(i, 1, n) rep(j, 1, m) a.num[j][i] = num[i][j]; return a; } //return the value of a determinate double val(){ int tot = 0; matrix tmp = (*this); if(n != m) { cout << "can't get the value of rectangle matrix, will return 0.0" << endl; return 0; } rep(i, 1, n) { if(abs(num[i][i]) < eps) { int cnt = i + 1; while(abs(num[cnt][i]) < eps && cnt <= n) cnt++; if(cnt == n + 1) return 0; else { swap_hang(i, cnt); tot++; } } rep(j, i + 1, n) { double cst = num[j][i] / num[i][i]; minus_lie(i, j, cst); } } double ans = 1; rep(i, 1, n) ans *= num[i][i]; *this = tmp; if(tot % 2) return -ans; return ans; } //return the submatrix of the matrix matrix submatrix(int ver1, int ver2, int hor1, int hor2){ matrix ans; ans.clear(); ans.n = hor2 - hor1 + 1; ans.m = ver2 - ver1 + 1; rep(i, 1, ans.n) rep(j, 1, ans.m) ans.num[i][j] = num[hor1 + i - 1][ver1 + j - 1]; return ans; } //if the solution is unique, return 1 and change the original n*(n+1) matrix into its solution //else return 0 bool la_solution1(){ if(abs(submatrix(1, n, 1, n).val()) < eps) return 0; rep(i, 1, n) { if(abs(num[i][i]) < eps) { int cnt = i + 1; while(abs(num[cnt][i]) < eps && cnt <= n) cnt++; swap_hang(i, cnt); } rep(j, i + 1, n) { double cst = num[j][i] / num[i][i]; minus_hang(i, j, cst); } } rep(i, 1, n) mul_hang(i, 1.0 / num[i][i]); repm(i, n, 2) { rep(j, 1, i - 1) minus_hang(i, j, num[j][i]); } return 1; } //convert the matrix into row echelon form matrix reduced(){ matrix x = (*this); int tot_cnt = 0; rep(i, 1, n) { if(i + tot_cnt > m) break; if(abs(x.num[i][i+tot_cnt]) < eps) { int cnt = -1; rep(j, i + 1, x.n) { if(abs(x.num[j][i+tot_cnt]) > eps) cnt = j; } if(cnt == -1) { tot_cnt++; i--; continue; } else { x.swap_hang(i, cnt); } } rep(j, i + 1, n) { double cst = x.num[j][i + tot_cnt] / x.num[i][i + tot_cnt]; x.minus_hang(i, j, cst); } } return x; } //convert the matrix into simplest form(only used for get the reverse) matrix normal_reduced() { matrix x = (*this).reduced(); rep(i, 1, n) { double tmp = x.num[i][i]; x.mul_hang(i, 1.0 / tmp); } repm(i, n, 2) { rep(j, 1, i - 1) x.minus_hang(i, j, x.num[j][i]); } return x; } //return the rank of the matrix int rank() { matrix x = (*this).reduced(); int i; for(i = 1; i <= x.n; i++) { int flag = 1; rep(j, 1, x.m) { if(abs(x.num[i][j]) > eps) flag = 0; } if(flag) break; } return i - 1; } //return the reverse of the matrix matrix reverse() { matrix x; if(abs(val()) < eps) { x = (*this); cout << "rank<(n-1),no reverse,will return the matrix" << endl; } x.clear(); x.n = n, x.m = 2 * n; rep(i, 1, n) rep(j, 1, n) x.num[i][j] = num[i][j]; rep(i, 1, n) x.num[i][n + i] = 1; x = x.normal_reduced(); x.print(); return x.submatrix(n + 1, 2 * n, 1, n); } //0 return the congruent standard form of the matrix B //1 return the matrix A, s.t. A.transposition*C*A==B (C is the original matrix) matrix normal_bi(bool flag) { matrix x; x.n = 2 * n; x.m = n; x.clear(); rep(i, 1, n) rep(j, 1, n) x.num[i][j] = num[i][j]; rep(i, 1, n) x.num[n + i][i] = 1; rep(i, 1, n - 1) { x.print(); if(abs(x.num[i][i]) < eps) { rep(j, i + 1, n) if(abs(x.num[j][i]) > eps) { x.minus_hang(j, i, -1); x.minus_lie(j, i, -1); break; } } if(abs(x.num[i][i]) < eps) { continue; } rep(j, i + 1, n) { double mul = x.num[j][i] / x.num[i][i]; x.minus_hang(i, j, mul); x.minus_lie(i, j, mul); } } if(flag) return x.submatrix(1, n, 1, n); else return x.submatrix(1, n, 1 + n, 2 * n); } /* //Jacobi method matrix diag() { matrix xx=(*this); double mx=0; do { int mi=0,mj=0; mx=0; rep(i,1,n) rep(j,1,n) { if(i!=j&&abs(num[i][j])>mx) { mx=abs(num[i][j]); mi=i;mj=j; } } if(mi>mj) swap(mi,mj); double ang; if(num[mi][mi]!=num[mj][mj]) ang=atan(-2.0*num[mi][mj]/(num[mi][mi]-num[mj][mj])); else { if(num[mi][mj]>0) ang=-pi/2; else ang=pi/2; } ang/=2.0; matrix tmp; tmp.m=tmp.n=n; tmp.normalize(); tmp.num[mi][mi]=tmp.num[mj][mj]=cos(ang); tmp.num[mj][mi]=sin(ang); tmp.num[mi][mj]=-sin(ang); (*this)=tmp.transposition()*(*this)*tmp; }while(mx>eps); swap(*this,xx); return xx; } */ }; //a class of matrix; the elements are polynomials struct poly_matrix{ int n, m; poly x[11][11]; poly_matrix() {} //generate a determinate, which value is the feature polynomial of the matrix orig poly_matrix(matrix orig): n(orig.n) { rep(i, 1, n) rep(j, 1, n) { vector<double> vec; vec.push_back(-orig.num[i][j]); if(i == j) vec.push_back(1.0); x[i][j] = poly(vec); } } //print the poly_matrix void print() { cout << endl; rep(i, 1, n) { rep(j, 1, n) x[i][j].print(); cout << endl; } cout << endl; } //swap the 2 lines of the poly_matrix void _swap(int ii, int jj) { rep(i, 1, n) swap(x[ii][i], x[jj][i]); } //update the jj(th) line with the value of line(jj)-fac*line(ii) void _sub(int ii, int jj, poly fac) { rep(i, 1, n) x[jj][i] = x[jj][i] - x[ii][i] * fac; } //update the ii(th) line with the value of line(ii)*fac void _mult(int ii, poly fac) { rep(i, 1, n) x[ii][i] = x[ii][i] * fac; } //return the value(a polynomial) of the poly_matrix poly feature_poly() { vector<double> tmp; tmp.push_back(1); poly_matrix ttmp = (*this); poly factor = poly(tmp); poly final = factor; rep(i, 1, n - 1) { int maxlen = 1000; int id = 0; rep(j, i, n) { if(!x[j][i].empty() && x[j][i].n < maxlen) { maxlen = x[j][i].n; id = j; } } if(id == 0) { vector<double> vec; vec.push_back(0.0); return poly(vec); } if(id != i) { factor = factor * -1.0; _swap(i, id); } rep(j, i + 1, n) { if(can_be_devided(x[j][i], x[i][i])) _sub(i, j, x[j][i] / x[i][i]); else { _mult(j, x[i][i]); factor = factor * x[i][i]; _sub(i, j, x[j][i] / x[i][i]); } } } rep(i, 1, n) final = final * x[i][i]; final = final / factor; (*this) = ttmp; return final; } }; //a class of vector struct vct{ int n; double num[101]; vct() {clear();} //set the ii(th) row as the vector's value vct(matrix x, int ii): n(x.n) { rep(i, 1, x.n) num[i] = x.num[i][ii]; clear(); } //return the dot of 2 vectors double friend operator * (vct a, vct b) { double ans = 0; rep(i, 1, a.n) ans += a.num[i] * b.num[i]; return ans; } //return the product of a vector and a number vct friend operator * (vct a, double b) { rep(i, 1, a.n) a.num[i] *= b; return a; } //return the addition of 2 vectors vct friend operator + (vct a, vct b) { rep(i, 1, a.n) a.num[i] += b.num[i]; return a; } //return the subtraction of 2 vectors vct friend operator - (vct a, vct b) { rep(i, 1, a.n) a.num[i] -= b.num[i]; return a; } //convert the vector into unitization form void normalize() { double sum = 0; rep(i, 1, n) sum += num[i] * num[i]; sum = sqrt(sum); rep(i, 1, n) num[i] /= sum; } //print the vector void print() { cout << endl; rep(i, 1, n) cout << num[i] << " "; cout << endl; } //clear the vector void clear() { rep(i, 0, 100) num[i] = 0; } }; //save the vector as the ii(th) row of the matrix void save(matrix &x, int ii, vct v) { rep(i, 1, v.n) x.num[i][ii] = v.num[i]; } bool cmp(double a,double b) { return abs(a) > abs(b); } //get the U, s.t. U^-1 * orig * U = diag, U is saved in as bool find_T(matrix x, matrix &as) { as.m = as.n = x.n; poly_matrix tmp = poly_matrix(x); poly p_tmp = tmp.feature_poly(); vector<double> vec = p_tmp.solu(); sort(vec.begin(), vec.end(), cmp); vector<double> feature; feature.push_back(vec[0]); rep(i, 1, vec.size() - 1) { if(abs(vec[i] - vec[i - 1]) > eps) feature.push_back(vec[i]); } vec = feature; matrix I; I.m = I.n = x.n; I.normalize(); vector<vct> vecc; rep(i, 0, vec.size() - 1) { matrix mat = I * vec[i] - x; mat = mat.reduced(); if(mat.rank() == 0) { rep(ii, 1, mat.n) { vct tmp; tmp.n = mat.n; tmp.num[ii] = 1.0; vecc.push_back(tmp); } continue; } bool main[101]; mem(main); rep(ii, 1, mat.rank()) { rep(j, 1, mat.m) { if(abs(mat.num[ii][j]) > eps) { main[j] = 1; break; } } } rep(ii, 1, mat.n) { if(!main[ii]) { matrix tmp; tmp.m = tmp.n = mat.rank(); int cnt = 0; rep(ti, 1, mat.n) { if(main[ti]) { cnt++; rep(j, 1, mat.rank()) tmp.num[j][cnt] = mat.num[j][ti]; } } tmp.m++; rep(j, 1, mat.rank()) tmp.num[j][tmp.m] = -mat.num[j][ii]; tmp.la_solution1(); vct ans; ans.n = mat.n; cnt = 0; rep(j, 1, mat.n) { if(main[j]) { cnt++; ans.num[j] = tmp.num[cnt][tmp.m]; } } ans.num[ii] = 1.0; vecc.push_back(ans); } } } if(vecc.size() < x.n) return 0; else { rep(i, 0, vecc.size() - 1) { vecc[i].normalize(); save(as, i + 1, vecc[i]); } } } /* void SVD(matrix x,matrix &s,matrix &v,matrix &d) { matrix tmp; tmp=x.transposition()*x; find_T(tmp,d); tmp=d.transposition()*x*d; int n=tmp.rank(); v=x; v.clear(); rep(i,1,n) v.num[i][i]=tmp.num[i][i]; } */
3281
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
