POJ 3526 The Teacher’s Side of Math
教你出题:给出根x = m√a + n√b,求原多项式方程。其中a和b是素数,m*n<20,多项式最高次数项的系数为1。
4.1更加复杂的数学问题
矩阵
首先题目的隐含条件是,多项式最高次数一定为m*n。将(a1/m + b1/n)k二项展开后,除了最高次数项被题目限定为1之外,各项的系数和必须为0。以各项系数为变量,列出线性方程组,然后高斯消元求解即可。
由于涉及到开方运算,所以使用了double和pow,涉及到精度问题,四舍五入糊弄过去了,就这么洒脱。
不过上次用的那个模板对题目给出的用例测试正常,但提交就WA。后来提交到了AOJ上,骗来了一份测试输入输出用例集,这才发现问题是上次那个算法精度损失太快了。于是换了个O(N3)的朴素算法,轻松AC~诸位可以自行下载The Teacher’s Side of Math 测试数据.zip。
#include <iostream> #include <algorithm> #include <vector> #include <cmath> using namespace std; #define LL long long #define MAX_N 20 + 8 int A, B, N, M; LL C[MAX_N][MAX_N]; // 生成组合数 void build_combinatorial_number() { for (int i = 0; i < MAX_N; ++i) { C[i][0] = C[i][i] = 1; for (int j = 1; j < i; ++j) C[i][j] = C[i - 1][j - 1] + C[i - 1][j]; } } inline int to_index(int n, int m) { return n * M + m + 1; } // 高斯消元法求解线性方程组 ax = b ,增广矩阵A=[a,b] vector<double> gauss(vector< vector<double> >& A) { int n = A.size(); for (int i = 0; i < n; ++i) { // Search for maximum in this column double maxEl = abs(A[i][i]); int maxRow = i; for (int k = i + 1; k < n; ++k) { if (abs(A[k][i]) > maxEl) { maxEl = abs(A[k][i]); maxRow = k; } } // Swap maximum row with current row (column by column) for (int k = i; k < n + 1; ++k) { double tmp = A[maxRow][k]; A[maxRow][k] = A[i][k]; A[i][k] = tmp; } // Make all rows below this one 0 in current column for (int k = i + 1; k < n; ++k) { double c = -A[k][i] / A[i][i]; for (int j = i; j < n + 1; ++j) { if (i == j) { A[k][j] = 0; } else { A[k][j] += c * A[i][j]; } } } } // Solve equation Ax=b for an upper triangular matrix A vector<double> x(n); for (int i = n - 1; i >= 0; --i) { x[i] = A[i][n] / A[i][i]; for (int k = i - 1; k >= 0; --k) { A[k][n] -= A[k][i] * x[i]; } } return x; } inline double round(double x) { return x >= 0.0 ? floor(x + 0.5) : ceil(x - 0.5); } void solve() { const int size = N * M + 1; vector< vector<double> > mat(size, vector<double>(size + 1, 0.0)); mat[0][0] = 1.0; mat[0][size] = 1.0; // the coefficient of its highest-order term (ad) is one,于是增广矩阵第一行肯定是[1, 0, 0,..., 1] for (int d = 0; d < size; ++d) { for (int k = 0; k <= d; ++k) { mat[to_index(k % N, (d - k) % M)][d < N * M ? d + 1 : 0] += C[d][k] * pow((double)A, k / N) * pow((double)B, (d - k) / M); } } vector<double> x = gauss(mat); reverse(x.begin() + 1, x.end()); for (int i = 0; i < size; ++i) { printf("%d%c", int(round(x[i])), i == N * M ? '\n' : ' '); } } ///////////////////////////SubMain////////////////////////////////// int main(int argc, char *argv[]) { #ifndef ONLINE_JUDGE freopen("in.txt", "r", stdin); freopen("out.txt", "w", stdout); #endif build_combinatorial_number(); while (~scanf("%d%d%d%d", &A, &N, &B, &M) && A > 0) { solve(); } #ifndef ONLINE_JUDGE fclose(stdin); fclose(stdout); system("out.txt"); #endif return 0; } ///////////////////////////End Sub//////////////////////////////////
14401901 | hankcs | 3526 | Accepted | 180K | 63MS | C++ | 2638B | 2015-07-15 22:32:45 |
Reference
http://martin-thoma.com/solving-linear-equations-with-gaussian-elimination/#tocAnchor-1-2
http://d.hatena.ne.jp/komiyam/20120924/1348415046
知识共享署名-非商业性使用-相同方式共享:码农场 » POJ 3526 The Teacher’s Side of Math 题解《挑战程序设计竞赛》
博主么么哒,您是acmer吗?
并不是