zoukankan      html  css  js  c++  java
  • 求矩阵的特征多项式

    一个比较慢的做法

      首先你要知道矩阵的特征多项式是什么。

      直接消元就可以了。

      时间复杂度:(O(n^5))(O(n^4))

    一个稍微快一点的做法

      观察到特征多项式的次数是(n)

      我们就可以插值。

      具体来说,先求出当(x=0ldots n)时特征多项式对应的点值,然后直接用拉格朗日插值插出来。

      时间复杂度:(O(n^4))

    一个更快的做法

      有一个性质:相似矩阵的特征多项式相同。

      所以可以把这个矩阵相似到一个可以快速求特征多项式的矩阵再求。

      怎么相似呢?

      (A)(B)相似意味着(A=PBP^{-1})

      假设当前矩阵是(A)。在消元过程中会左乘一个初等矩阵(P),这时在(A)的右边乘上这个初等矩阵的逆矩阵(P^{-1}),得到(PAP^{-1})

      显然这个矩阵和矩阵(A)相似。

      但是我们只能把所有(igeq j+2)的部分消成(0),也就是说会消成一个上海森堡矩阵。

      接下来就要快速求一个上海森堡矩阵的特征多项式。

      上海森堡矩阵长这样:

    [egin{pmatrix} a_{1,1}&a_{1,2}&a_{1,3}&cdots&a_{1,n}\ a_{2,1}&a_{2,2}&a_{2,3}&cdots&a_{2,n}\ 0&a_{3,2}&a_{3,3}&cdots&a_{3,n}\ vdots&vdots&vdots&ddots&vdots\ 0&0&0&cdots&a_{n,n} end{pmatrix} ]

      设右下角(i imes i)的矩阵的特征多项式为(f_i)

      考虑对第一列展开。

      那么肯定是连续的消掉若干个第二行后再消掉第一行。

      即

    [f_i=a_{i,i}f_{i+1}-a_{i+1,i}a_{i,i+1}f_{i+2}+a_{i+1,i}a_{i+2,i+1}a_{i,i+2}f_{i+3}+cdots ]

      注意到只有主对角线上的元素一次的,其他都是常数的,所以前面的系数的次数最多是(1)

      这样可以从后往前求出所有(f_i)。时间复杂度是(O(n^3))

      这样我们就得到了一个时间复杂度是(O(n^3))的优秀做法。

      这个做法有什么用?

      可以去卡别人。比如用(O(n^3+n^2log m))的做法卡(O(n^3log m))的矩阵快速幂。

    代码

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<utility>
    using namespace std;
    typedef long long ll;
    typedef pair<int,int> pii;
    typedef pair<ll,ll> pll;
    const ll p=998244353;
    struct poly
    {
        ll a[510];
        int n;
        poly()
        {
            memset(a,0,sizeof a);
            n=0;
        }
        ll &operator [](int x)
        {
            return a[x];
        }
    };
    int n;
    ll a[510][510];
    void add(poly &a,poly &b,pll c)
    {
        a.n=max(a.n,b.n+bool(c.first));
        int i;
        for(i=0;i<=a.n;i++)
        {
            a[i]=(a[i]+b[i]*c.second)%p;
            if(i)
                a[i]=(a[i]+b[i-1]*c.first)%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;
    }
    void gao1()
    {
        int i,j,k;
        for(i=1;i<=n;i++)
        {
            for(j=i+1;j<=n;j++)
                if(a[i][j])
                    break;
            if(j>n)
                continue;
            if(j!=i+1)
            {
                for(k=i;k<=n;k++)
                    swap(a[i+1][k],a[j][k]);
                for(k=1;k<=n;k++)
                    swap(a[k][i+1],a[k][j]);
            }
            ll e=fp(a[i+1][i],p-2);
            for(j=i+2;j<=n;j++)
                if(a[j][i])
                {
                    ll v=e*a[j][i]%p;
                    for(k=i;k<=n;k++)
                        a[j][k]=(a[j][k]-a[i+1][k]*v)%p;
                    for(k=1;k<=n;k++)
                        a[k][i+1]=(a[k][i+1]+a[k][j]*v)%p;
                }
        }
    }
    pll c[510][510];
    poly f[510];
    pll operator *(pll a,pll b)
    {
        return pll((a.first*b.second+a.second*b.first)%p,a.second*b.second%p);
    }
    pll operator *(pll a,ll b)
    {
        return pll(a.first*b%p,a.second*b%p);
    }
    void gao2()
    {
        f[n+1][0]=1;
        int i,j;
        for(i=n;i>=1;i--)
        {
            pll v(0,1);
            for(j=i+1;j<=n+1;j++)
            {
                add(f[i],f[j],v*c[i][j-1]*((j-i)&1?1:-1));
                v=v*c[j][j-1];
            }
        }
    }
    int main()
    {
        scanf("%d",&n);
        int i,j;
        for(i=1;i<=n;i++)
            for(j=1;j<=n;j++)
                scanf("%lld",&a[i][j]);
        gao1();
        for(i=1;i<=n;i++)
            for(j=1;j<=n;j++)
            {
                c[i][j].second=-a[i][j];
                if(i==j)
                    c[i][j].first=1;
            }
        gao2();
    	for(i=0;i<=n;i++)
    		printf("%lld ",(f[1][i]+p)%p);
    	printf("
    ");
        return 0;
    }
    
  • 相关阅读:
    表的相关操作
    存储引擎介绍
    库的相关操作
    初始数据库
    linux版本的mysql安装
    mysql在windows上的安装即配置
    线程实际操作篇
    用户模板和用户场景
    顶会热词统计
    移动端疫情展示
  • 原文地址:https://www.cnblogs.com/ywwyww/p/8522541.html
Copyright © 2011-2022 走看看