放牧代码和思想
专注自然语言处理、机器学习算法
    This thing called love. Know I would've. Thrown it all away. Wouldn't hesitate.

POJ 3526 The Teacher’s Side of Math 题解《挑战程序设计竞赛》

目录

POJ 3526 The Teacher’s Side of Math 

教你出题:给出根x = ma + nb,求原多项式方程。其中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 题解《挑战程序设计竞赛》

评论 2

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址
  1. #1

    博主么么哒,您是acmer吗?

    逍遥小章9年前 (2015-07-27)回复

我的作品

HanLP自然语言处理包《自然语言处理入门》