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~varint 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 cinline 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;}