Matrix Power Series(POJ3233 矩阵快速幂+二分)

链接:link.
题意:给一个n * n的矩阵A,求S=A+A2+A3+…+Ak;
思路:构造
( E A O A ) \begin{pmatrix}{E}&{A}\\ {O}&{A}\\\end{pmatrix} (EOAA) * ( S k − 1 A k − 1 ) \begin{pmatrix}{S_{k-1}}\\ {A^{k-1}}\\\end{pmatrix} (Sk1Ak1)= ( S k A k ) \begin{pmatrix}{S_{k}}\\ {A^{k}}\\\end{pmatrix} (SkAk)
二分:A+A2+A3+A4=A2 * (A+A2);

int n,k,mod;
struct node
{
    int v[35][35];
}res,init;
node mult(node x,node y)
{
    node tmp;
    memset(tmp.v,0,sizeof(tmp.v));
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            for(int k=0;k<n;k++)
            {
                tmp.v[i][j]=(tmp.v[i][j]+x.v[i][k]*y.v[k][j]%mod)%mod;
            }
        }
    }
    return tmp;
}
node Pow(node x,int k)
{
    node tmp;
    memset(tmp.v,0,sizeof(tmp.v));
    for(int i=0;i<n;i++)tmp.v[i][i]=1;
    while(k)
    {
        if(k&1)tmp=mult(tmp,x);
        x=mult(x,x);
        k=k>>1;
    }
    return tmp;
}
node add(node x,node y)
{
    node tmp;
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            tmp.v[i][j]=(x.v[i][j]%mod+y.v[i][j]%mod)%mod;
        }
    }
    return tmp;
}
node sum(node x,int k)
{
    if(k==1)return x;
    node tmp=sum(x,k/2);
    if(k&1)
    {
        node y=Pow(x,k/2+1);
        tmp=add(mult(y,tmp),tmp);
        return add(tmp,y);
    }
    else
    {
        node y=Pow(x,k/2);
        return add(mult(y,tmp),tmp);
    }
}
int main()
{
    cin>>n>>k>>mod;
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            int x;cin>>x;
            init.v[i][j]=x%mod;
        }
    }
    res=sum(init,k);
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            cout<<res.v[i][j]%mod<<" ";
        }
        cout<<"\n";
    }
    return 0;
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值