【poj3233】Matrix Power Series（递推+矩阵快速幂）

最新推荐文章于 2022-05-05 20:34:22 发布

原创最新推荐文章于 2022-05-05 20:34:22 发布 · 855 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#DP大法 #矩阵快速幂 #poj

题解同时被 3 个专栏收录

1097 篇文章

订阅专栏

245 篇文章

订阅专栏

矩阵

21 篇文章

订阅专栏

该博客介绍了如何解决POJ3233问题，重点讨论了利用递推关系和矩阵快速幂优化算法以提高效率，而非采用朴素方法。博主分享了他们的思考过程和代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目描述

传送门

题解

用朴素的矩阵快速幂的话时间无法承受。其实这道题的思路挺奇特的，不是自己想出来的很不爽。
这里有一个递推的思路：考虑k能不能从之前的状态转移过来。答案是肯定的。
当 $k\%2==0$ 时：

s u m (k) = \sum i = 1 k A k = (1 + A k 2) * (A + A 2 + \dots + A k 2)

$sum(k)=\sum_{i=1}^k A^k=(1+A^{k\over 2})*(A+A^2+…+A^{k\over 2})$
当

k%2!=0 $k\%2!=0$ 时

s u m (k) = \sum i = 1 k A k = A + (A + A ⌊ k 2 ⌋) * (A + A 2 + \dots + A ⌊ k 2 ⌋)

$sum(k)=\sum_{i=1}^k A^k=A+(A+A^{\lfloor {k\over 2}\rfloor})*(A+A^2+…+A^{\lfloor {k\over 2}\rfloor})$
也就是说先计算出

sum(⌊k2⌋) $sum(\lfloor {k\over 2}\rfloor)$ 的值就可以经过很少的计算得出

sum(k) $sum(k)$ 的值，那么我们可以不断递归下去。
有必要说明一点：第一个公式出现了1，它其实就代表了一个单位矩阵，即

I = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 100 ⋮ 0 010 ⋮ 0 001 ⋮ 0 \dots \dots \dots ⋱ \dots 000 ⋮ 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$I= \begin{bmatrix} 1 & 0 & 0&\cdots&0 \\ 0 & 1 & 0&\cdots&0\\ 0 & 0 & 1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&1\\ \end{bmatrix}$
满足

A∗I=I∗A=A $A*I=I*A=A$ （

A 为任意矩阵）

代码

#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;

struct hp{
    int a[35][35];
}A,unit;
int n,k,m;

inline hp jia(hp a,hp b){
    hp c;
    for (int i=1;i<=n;++i)
      for (int j=1;j<=n;++j)
        c.a[i][j]=(a.a[i][j]+b.a[i][j])%m;
    return c;
}

inline hp cheng(hp a,hp b){
    hp c;
    memset(c.a,0,sizeof(c.a));
    for (int i=1;i<=n;++i)
      for (int j=1;j<=n;++j)
        for (int k=1;k<=n;++k)
          c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%m)%m;
    return c;
}

inline hp fast_pow(hp a,int p){
    hp ans=unit;
    for (;p;p>>=1,a=cheng(a,a))
      if (p&1)
        ans=cheng(ans,a);
    return ans;
}

inline hp doit(int k){
    if (k==1) return A;
    hp now=doit(k/2);
    if (k%2==0){
        hp x=fast_pow(A,k/2);
        x=jia(unit,x);
        x=cheng(x,now);
        return x;
    }
    else{
        hp x=fast_pow(A,k/2+1);
        x=jia(A,x);
        x=cheng(x,now);
        x=jia(A,x);
        return x;
    }
}

int main(){
    scanf("%d%d%d",&n,&k,&m);
    for (int i=1;i<=n;++i)
      for (int j=1;j<=n;++j)
        scanf("%d",&A.a[i][j]);
    for (int i=1;i<=n;++i)
      unit.a[i][i]=1;
    hp ans=doit(k);
    for (int i=1;i<=n;++i)
      for (int j=1;j<=n;++j)
        printf("%d%c",ans.a[i][j]," \n"[j==n]);
}