【XSY3154】入门多项式 高斯消元

题目大意

　　给你一个 \(n\times n\)的矩阵 \(A\)，求次数最小且最高次项为 \(1\) 的多项式 \(F(x)\)，满足 \(F(A)=0\)。

　　所有操作都对 \(p\) 取模。

　　\(n\leq 70,n<p\leq 998244353\)

题解

　　显然特征多项式满足条件，但不一定是最优的。

　　设答案为 \(F(x)=\sum_{i\geq 0}f_ix^i\)。

　　那么
\[ \begin{cases} f_0{(A^0)}_{1,1}+f_1{(A^1)}_{1,1}+\cdots+f_n{(A^n)}_{1,1}&=0\\ f_0{(A^0)}_{1,2}+f_1{(A^1)}_{1,2}+\cdots+f_n{(A^n)}_{1,2}&=0\\ \vdots\\ f_0{(A^0)}_{n,n}+f_1{(A^1)}_{n,n}+\cdots+f_n{(A^n)}_{n,n}&=0 \end{cases} \]
　　这就是一个方程组，可以通过高斯消元来求解。

　　观察高斯消元的过程。

　　如果在消第 \(i\) 列的时候找不到主元，就说明这个矩阵的前 \(i\) 列不满秩，那么就可以钦定 \(f_{i-1}=1\)，从而得到一组解。

　　否则前 \(i\) 列是满秩的，唯一可能的解为 \(f_0=f_1=\ldots=f_{i-1}=0\)

　　时间复杂度：\(O(n^4)\)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=80;
int n;
ll p;
ll fp(ll a,ll b)
{
    ll s=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)
            s=s*a%p;
    return s;
}
struct mat
{
    ll a[N][N];
    mat()
    {
        memset(a,0,sizeof a);
    }
    ll *operator [](int x)
    {
        return a[x];
    }
};
mat operator *(mat a,mat b)
{
    mat c;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
        {
            __int128 s=0;
            for(int k=1;k<=n;k++)
                s+=(ll)a[i][k]*b[k][j];
            c[i][j]=s%p;
        }
    return c;
}
mat a[N];
ll ans[N];
ll c[N*N][N];
int m;
void gao(int x)
{
    ans[x]=1;
    for(int i=1;i<x;i++)
        ans[i]=(-c[i][x]*fp(c[i][i],p-2)%p+p)%p;
    printf("%d\n",x-1);
    for(int i=1;i<=x;i++)
        printf("%lld ",ans[i]);
}
void gao()
{
    for(int i=1;i<=n+1;i++)
    {
        int flag=0;
        for(int j=i;j<=m;j++)
            if(c[j][i])
            {
                flag=j;
                break;
            }
        if(!flag)
        {
            gao(i);
            return;
        }
        if(flag!=i)
        {
            for(int k=i;k<=n+1;k++)
                swap(c[i][k],c[flag][k]);
        }
        ll inv=fp(c[i][i],p-2);
        for(int j=1;j<=m;j++)
            if(j!=i&&c[j][i])
            {
                ll v=c[j][i]*inv%p;
                for(int k=i;k<=n+1;k++)
                    c[j][k]=(c[j][k]-v*c[i][k])%p;
            }
    }
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("a.in","r",stdin);
    freopen("a.out","w",stdout);
#endif
    scanf("%d%lld",&n,&p);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            scanf("%lld",&a[1][i][j]);
    for(int i=1;i<=n;i++)
        a[0][i][i]=1;
    for(int i=2;i<=n;i++)
        a[i]=a[i-1]*a[1];
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
        {
            m++;
            for(int k=0;k<=n;k++)
                c[m][k+1]=a[k][i][j];
        }
    gao();
    return 0;
}

转载于:https://www.cnblogs.com/ywwyww/p/9279448.html