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

Posted 2020-10-26 Styx-ferryman

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3
题意：给你一个矩阵A，且已知S = A + A^2 + A^3 +...+ A^k，求S
题解：假设有A + A^2 + A^3 +...+ A^n
那么如果n为偶数，该式子可以写为
（I+A^(N/2)）*(A+A^2+...+A^(N/2))
如果n为奇数，则可以拆成一个矩阵A^n和n为偶数的另一串式子。
于是就可以分治了，1700ms莫名慌
学到了如何memset结构体里的数组，还是有收获的。
代码如下：

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

struct matrix
{
    int a[35][35];
    matrix()
    {
        memset(a,0,sizeof(a));
    }
};

int n,k,m;

void print(matrix &ans)
{
    for(int i=1; i<=n; i++)
    {
        for(int j=1; j<=n-1; j++)
        {
            printf("%d ",ans.a[i][j]);
        }
        printf("%d\n",ans.a[i][n]);
    }
}

matrix add(matrix a,matrix b)
{
    matrix 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];
            if(c.a[i][j]>=m)
            {
                c.a[i][j]%=m;
            }
        }
    }
    return c;
}

matrix mul(matrix a,matrix b)
{
    matrix c;
    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]+=a.a[i][k]*b.a[k][j];
                if(c.a[i][j]>=m)
                {
                    c.a[i][j]%=m;
                }
            }
        }
    }
    return c;
}

matrix kasumi(matrix a,int b)
{
    matrix ans;
    for(int i=1; i<=n; i++)
    {
        ans.a[i][i]=1;
    }
    if(b==0)
    {
        return ans;
    }
    if(b==1)
    {
        return a;
    }
    while(b)
    {
        if(b&1)
        {
            ans=mul(ans,a);
        }
        a=mul(a,a);
        b>>=1;
    }
    return ans;
}

matrix solve(matrix a,int k)
{
    if(k==1)
    {
        return a;
    }
    if(k&1)
    {
        return add(solve(a,k-1),kasumi(a,k));
    }
    else
    {
        return mul(add(kasumi(a,0),kasumi(a,k>>1)),solve(a,k>>1));
    }
}

int main()
{
    scanf("%d%d%d",&n,&k,&m);
    matrix a,ans;
    for(int i=1; i<=n; i++)
    {
        for(int j=1; j<=n; j++)
        {
            scanf("%d",&a.a[i][j]);
        }
    }
    ans=solve(a,k);
    print(ans);
}