Matrix Power Series(矩阵加速，矩阵套矩阵)

最新推荐文章于 2023-08-16 10:14:07 发布

霖_麒

最新推荐文章于 2023-08-16 10:14:07 发布

阅读量705

点赞数

分类专栏：数论题解

本文链接：https://blog.csdn.net/zyz_bz/article/details/89311150

版权

题解同时被 2 个专栏收录

25 篇文章 0 订阅

订阅专栏

数论

14 篇文章 1 订阅

订阅专栏

文章目录

不知道矩阵的点这里

题目

描述
给定矩阵 $A$ ，求矩阵 $S=A^1+A^2+…+A^k$ ，输出矩阵， $S$ 矩阵中每个元都要模 $m$ 。
数据范围： $n (n ≤ 30) , k (k ≤ 10^9) ，m (m < 10^4)$
输入
输入三个正整数 $n ， k ， m$
输出
输出矩阵 $\%m$
样例输入
2 2 4
0 1
1 1
样例输出
1 2
2 3

思路

我们从一般到特殊，我们先想，若 $a$ 仅仅是一个数，我们如何计算 $a^1+a^2+…+a^k$ ，我相信肯定有人喊：等比数列。但请注意：我们要用矩阵求解
第一，我们肯定要拿个值存 $a^1+a^2+…+a^k$ ，所以我们定 $S_k=a^1+a^2+……+a^k$ $S_{k+1}=a^1+a^2+…+a^k+a^{k+1}$ $S_k=a^1+a^2+……+a^k$ $S_{k+1}=S_k*a+a（看的出来不？）=S_k+a^{k+1}$ 所以根据第一个等式就有了下面矩阵
$\begin{bmatrix} S_k ，1 \end{bmatrix}* \begin{bmatrix} a ，0\\ a，1 \end{bmatrix}$
而继续根据第二个等式就可以再得出一个矩阵
$\begin{bmatrix} S_k ，a^k \end{bmatrix}* \begin{bmatrix} 1 ，0\\ a，a \end{bmatrix}$
如何优化成一个矩阵的次幂而不去乘另一个矩阵呢
首先 $\begin{bmatrix} a^k ，S_k\end{bmatrix}$ 这个矩阵肯定是要的，我们要变成 $\begin{bmatrix} a^{k+1} ，S_{k+1}\end{bmatrix}$
因为矩阵是第一行乘第一列，所以 $a^k$ 要乘 $a$ ，而 $S_k$ 不存在了，所以 $S_k$ 要乘0
而第一行乘第二列得到 $S_{k+1}$ ，因为 $S_{k+1}=S_k+a^{k+1}$ ，所以 $a^k$ 要乘 $a$ ，而 $S_k$ 也要加上，所以 $S_k$ 要乘1，所以就很清楚了
$\begin{bmatrix} a^k ，S_{k-1}\end{bmatrix}* \begin{bmatrix} a ，1(S_0)\\ 0，1 \end{bmatrix}= \begin{bmatrix} a^{k+1} ，S_{k}\end{bmatrix}$ 而我们观察到第一个矩阵的第一行和第二个矩阵的第一行的形式是一样的，所以我们可以合并成 ${\begin{bmatrix} a ，1\\0，1\end{bmatrix}}^{k+1}$ 它的第一行的第二列就是答案了
把 $a$ 推广成矩阵 $A$ 就行了
这两个矩阵请读者自行尝试（注意细节）
$\begin{bmatrix} S_k ，1 \end{bmatrix}* \begin{bmatrix} a ，0\\ a，1 \end{bmatrix}$

$\begin{bmatrix} S_k ，a^k \end{bmatrix}* \begin{bmatrix} 1 ，0\\ a，a \end{bmatrix}$

代码

我用的第三个矩阵
再提醒一点，注意细节

#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define LL long long
 
int n, m, k;
 
struct Mx{
    int x, y;
    LL a[33][33];
 
    inline void cl(){
        memset(a, 0, sizeof(a));
    }
 
    Mx operator * (const Mx &b)const{
        Mx r;
        for(int i = 1; i <= x; i ++){
            for(int j = 1; j <= b.y; j ++){
                r.a[i][j] = 0;
                for(int h = 1; h <= b.x; h ++)
                    r.a[i][j] = (r.a[i][j]+b.a[h][j]*a[i][h]%m)%m;
            }
        }
        r.x = x;
        r.y = b.y;
        return r;
    }
 
    Mx operator + (const Mx &b)const{
        Mx r;
        r.cl();
        for(int i = 1; i <= x; i ++){
            for(int j = 1; j <= y; j ++)
                r.a[i][j] = (r.a[i][j]+b.a[i][j]+a[i][j])%m;
        }
        r.x = b.x;
        r.y = b.y;
        return r;
    }
 
    Mx operator - (const Mx &b)const{
        Mx r;
        r.cl();
        for(int i = 1; i <= x; i ++){
            for(int j = 1; j <= y; j ++)
                r.a[i][j] = (r.a[i][j]-b.a[i][j]+a[i][j]+m*m)%m;
        }
        r.x = b.x;
        r.y = b.y;
        return r;
    }
 
    void print(){
        for(int i = 1; i <= x; i ++){
            for(int j = 1; j < y; j ++)
                printf("%lld ",a[i][j]);
            printf("%lld\n", a[i][y]);
        }
    }
}one;
 
struct Mxx{
    int xx, yy;
    Mx A[3][3];
 
    Mxx operator * (const Mxx &b)const{
        Mxx r;
        for(int i = 1; i <= xx; i ++){
            for(int j = 1; j <= b.yy; j ++){
                r.A[i][j].cl();
                r.A[i][j].x = r.A[i][j].y = n;
                for(int h = 1; h <= b.xx; h ++)
                    r.A[i][j] = r.A[i][j] + (A[i][h] * b.A[h][j]);
            }
        }
        r.xx = xx;
        r.yy = b.yy;
        return r;
    }
 
}F;
 
inline Mxx qkp(Mxx x, int y){
    Mxx sum;
    sum.xx = sum.yy = x.xx;
    for(int i = 1; i <= x.xx; i ++)
        sum.A[i][i] = one;
    sum.A[2][1].cl();
    sum.A[1][2].cl();
    sum.A[1][1].x = sum.A[1][1].y = n;
    sum.A[1][2].x = sum.A[1][2].y = n;
    sum.A[2][2].x = sum.A[2][2].y = n;
    sum.A[2][1].x = sum.A[2][1].y = n;
    while( y ){
        if( y&1 ){
            sum = x * sum;
        }
        x = x * x;
        y >>= 1;
    }
    return sum;
}
 
int main(){
    scanf("%d%d%d",&n,&k,&m);
    F.xx = F.yy = 2;
    F.A[1][1].cl();
    F.A[2][1].cl();
    F.A[1][2].cl();
    F.A[2][2].cl();
    F.A[1][1].x = F.A[1][1].y = n;
    F.A[1][2].x = F.A[1][2].y = n;
    F.A[2][2].x = F.A[2][2].y = n;
    F.A[2][1].x = F.A[2][1].y = n;
    one.cl();
    one.x = one.y = n;
    for(int i = 1; i <= n; i ++){
        one.a[i][i] = 1;
        for(int j = 1; j <= n; j ++)
            scanf("%lld",&F.A[1][1].a[i][j]);
    }
    F.A[1][2] = F.A[2][2] = one;
    F = qkp(F, k+1);
    F.A[1][2] = F.A[1][2] - one;
    F.A[1][2].print();
    return 0;
}