POJ2065 SETI 高斯消元

题目链接：http://poj.org/problem?id=2065

题目大意：给你一个素数P(P<=30000)和一串长为n的字符串str[]。字母'*'代表0，字母a-z分别代表1-26，这n个字符所代表的数字分别代表f(1)、f(2)....f(n)。

定义: f (k) = ∑0<=i<=n-1aiki (mod p) (1<=k<=n,0<=ai<P)

求a0、a1.....an-1。题目保证肯定有唯一解。

分析：高斯消元。根据上面的公式显然可以列出有n个未知数的n个方程式：

a0*1^0 + a1*1^1+a2*1^2+........+an-1*1^(n-1) = f(1)

a0*2^0 + a1*2^1+a2*2^2+........+an-1*2^(n-1) = f(2)

......

a0*n^0 + a1*n^1+a2*n^2+........+an-1*n^(n-1) = f(n)

然后采用高斯消元法来解上面的方程组即可。

实现代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
const int maxn=70;
int mod,a[maxn][maxn+1],x[maxn],equ,var;
char s[100];
void init()
{
    equ=strlen(s);
    var=equ;
    for(int i=0;i<var;i++)
    {
        int tmp=1;
        for(int j=0;j<var;j++)
        {
           a[i][j]=tmp;
           tmp=tmp*(i+1)%mod;
        }
        a[i][var]= (s[i]=='*')?0:(s[i]-'a'+1);
    }
}
void debug()
{
    for(int i=0;i<equ;i++)
    {
        for(int j=0;j<=var;j++)
          cout<<a[i][j]<<" ";
        cout<<endl;
    }
    cout<<endl;
}
int gcd(int a,int b)
{
    if(a<0) return gcd(-a,b);
    if(b<0) return gcd(a,-b);
    return b==0?a:gcd(b,a%b);
}
int ex_gcd(int a,int b,int &x,int &y)
{
    if(b==0){
         x=1,y=0;
         return a;
    }
    int ret=ex_gcd(b,a%b,x,y);
    int tmp=x;
    x=y;
    y=tmp-a/b*y;
    return ret;
}
void gauss()
{
    int k,col;
    for(k=col=0;k<equ&&col<var;k++,col++)
    {
        int mx=k;
        for(int i=k+1;i<equ;i++)
           if(a[i][col]>a[mx][col]) mx=col;
        if(k!=col)
        {
            for(int i=k;i<var+1;i++)
               swap(a[k][i],a[mx][i]);
        }
        if(!a[k][col])
        {
            k--;continue;
        }
        for(int i=k+1;i<equ;i++)
        if(a[i][col])
        {
            int lcm=a[k][col]/gcd(a[k][col],a[i][col])*a[i][col];
            int ta=lcm/a[i][col],tb=lcm/a[k][col];
            for(int j=col;j<var+1;j++)
               a[i][j]=((a[i][j]*ta-a[k][j]*tb)%mod+mod)%mod;
        }
    }
    for(int i=k;i<equ;i++)
      if(a[i][col]) return ;

    for(int i=0,j;i<equ;i++)
       if(!a[i][i])
       {
           for(j=i+1;j<var;j++)
             if(a[i][j]) break;
           if(j>=var) break;
           for(int r=0;r<equ;r++) swap(a[r][i],a[r][j]);
       }

    // debug();
    for(int i=k-1;i>=0;i--)
    {
        int tmp=a[i][var],x1,y1;
        for(int j=i+1;j<var;j++)
          tmp=((tmp-a[i][j]*x[j])%mod+mod)%mod;
        int gcd=ex_gcd(a[i][i],mod,x1,y1);
        x[i]=(x1*tmp/gcd%mod+mod)%mod;
    }
    printf("%d",x[0]);
    for(int i=1;i<var;i++) printf(" %d",x[i]);
    puts("");
}
int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d %s",&mod,s);
        init();
      //  debug();
        gauss();
    }
    return 0;
}