poj 2065 SETI（高斯消元）

http://poj.org/problem?id=2065

题意：

输入一个素数p和一个字符串s（只包含小写字母和‘*’）,字符串中每个字符对应一个数字，'*'对应0，‘a’对应1，‘b’对应2....

例如str[] = "abc", 那么说明 n=3, 字符串所对应的数列为1, 2, 3。

题目中定义了一个函数：

a0*1^0 + a1*1^1+a2*1^2+........+an-1*1^(n-1) = f(1)(mod p)， f(1) = str[0] = a = 1;
a0*2^0 + a1*2^1+a2*2^2+........+an-1*2^(n-1) = f(2)(mod p)， f(2) = str[1] = b = 2;
..........
a0*n^0 + a1*n^1+a2*n^2+........+an-1*n^(n-1) = f(n)(mod p)，f(n) = str[n-1] = ````

求出 a0,a1,a2....an-1.

思路:若字符串长度为n，那么这是一个含有n个方程n个未知数的线性方程组。所以只需把相应的系数转为成矩阵解方程组。高斯消元+同余方程。。

刚学的高斯消元，之前学的线性代数也忘的差不多了，又回头瞄了瞄。

#include <queue>
#include <algorithm>
#include <stack>
#include <vector>
#include <string.h>
#include <stdio.h>
using namespace std;

int p;
char s[80];
int equ,var; //有equ个方程，var个变元。增光矩阵行数是equ，从0~equ-1，列数是var+1,从0~var
int a[80][80];
int x[80];//解集

int gcd(int a, int b)
{
	if(b == 0)
		return a;
	return gcd(b,a%b);
}

inline int lcm(int a, int b)
{
	return a*b / gcd(a,b);
}

//快速幂取模 a^b mod c
inline int mod_exp(int a, int b, int c)
{
	int res = 1;
	int t = a%c;

	while(b)
	{
		if(b&1)
			res = res*t%c;
		t = t*t%c;
		b >>= 1;
	}
	return res;
}

int Gauss()
{
	int i,j,row,col,max_r;
	int ta,tb,tmp,l;

	row = col = 0;

	while(row < equ && col < var)
	{
		max_r = row; //在col列，max_r初始化为当前行最大
		for(i = row+1; i < equ; i++)
		{
			if( abs(a[i][col]) > abs( a[max_r][col]) )
				max_r = i;
		}
		//找到最大行与当前行交换
		if(max_r != row)
		{
			for(j = col; j < var+1; j++)
				swap(a[row][j], a[max_r][j]);
		}
		//交换后如果该行最大值是0，直接处理该行的下一列
		if(a[row][col] == 0)
		{
			col++;
			continue;
		}
		
		for(i = row+1; i < equ; i++)
		{
			if(a[i][col] == 0) continue;
			l = lcm( abs(a[i][col]), abs(a[row][col]) );
			ta = l/a[i][col];
			tb = l/a[row][col];
			if(a[i][col] * a[row][col] < 0) tb = -tb; //注意异号

			for(j = col; j < var+1; j++)
				a[i][j] = ( (a[i][j]*ta - a[row][j]*tb)%p + p)%p;
		}
		row++;
		col++;
	}
	// 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).  
	for(i = row; i < var; i++)
		if(a[i][var] != 0)
			return -1;
	// 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.返回自由变元的个数
	if(row < var)
		return var-row;
	// 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵. 
	// 计算出Xn-1, Xn-2 ... X0.   
	for(i = var-1; i >= 0; i--)
	{
		tmp = a[i][var];
		for(j = i+1; j < var; j++)
			tmp = ( (tmp-a[i][j]*x[j]) %p + p) %p;
		while( tmp % a[i][i] != 0)
			tmp += p;
		x[i] = ( tmp / a[i][i] ) % p;
	}
	return 0;
}

int main()
{
	int test;
	scanf("%d",&test);
	while(test--)
	{
		scanf("%d %s",&p,s);

		equ = var = strlen(s);

		for(int i = 0; i < equ; i++)
		{
			if(s[i] == '*')
				a[i][var] = 0;
			else a[i][var] = s[i]-'a'+1;

			for(int j = 0; j < var; j++)
				a[i][j] = mod_exp(i+1,j,p);
		}
	
		Gauss();

		for(int i = 0; i < var-1; i++)
			printf("%d ",x[i]);
		printf("%d\n",x[var-1]);

	}
	return 0;
}