POJ 3233 Matrix Power Series (矩阵快速幂 + 二分思想）

最新推荐文章于 2019-09-07 19:49:23 发布

原创最新推荐文章于 2019-09-07 19:49:23 发布 · 392 阅读

0 ·

CC 4.0 BY-SA版权

本文介绍了一种利用矩阵快速幂高效计算多项式求和的方法。通过对矩阵进行特定的组合和快速幂运算，可以在log级别的时间复杂度内完成计算。适用于解决数据范围大的问题。

题目连接

题意：求 A + A² + A³ + … + A^k.

数据范围：（n (n ≤ 30), k (k ≤ 10⁹) and m (mod< 10⁴)

根据矩阵的分配结合率可以把上式如下处理：

假设k=6，A^1+A^2+A^3+A^4+A^5+A^6 = A^1+A^2+A^3 + A^3*(A^1+A^2+A^3) = (E + A^3)*(A^1+A^2+A^3);

而，A^1+A^2+A^3 = (A^1+E)*A^1 + A^3;

所以 A^1+A^2+A^3+A^4+A^5+A^6 = （A^3+E ) * [ (A^1+E)*A^1 + A^3 ]

可以看到，每次k除2，最多只会有logk个形如（A^i + E）这样的式子相乘，然后如果k不能整除2的话，就把最后一个拿出来单独加，整个复杂度是logk的，求A^k次是n^3*logk,总复杂度n^3*logk*logk。

【代码】

/* ***********************************************
Author        :angon

************************************************ */
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <stack>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;
#define showtime fprintf(stderr,"time = %.15f\n",clock() / (double)CLOCKS_PER_SEC)
#define lld %I64d
#define REP(i,k,n) for(int i=k;i<n;i++)
#define REPP(i,k,n) for(int i=k;i<=n;i++)
#define scan(d) scanf("%d",&d)
#define scanl(d) scanf("%I64d",&d)
#define scann(n,m) scanf("%d%d",&n,&m)
#define scannl(n,m) scanf("%I64d%I64d",&n,&m)
#define mst(a,k)  memset(a,k,sizeof(a))
#define LL long long
#define N 30
inline int read(){int s=0;char ch=getchar();for(; ch<'0'||ch>'9'; ch=getchar());for(; ch>='0'&&ch<='9'; ch=getchar())s=s*10+ch-'0';return s;}

int n,k,mod;

struct Matrix
{
    int v[N][N];
};
Matrix operator * (Matrix A, Matrix B)
{
    Matrix ans;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            ans.v[i][j] = 0;
            for (int k = 0; k < n; k ++)
            {
                ans.v[i][j] += (A.v[i][k] * B.v[k][j]) % mod;
            }
            ans.v[i][j] %= mod;
        }
    }
    return ans;
}

Matrix matrix_pow(Matrix C, int p)
{
    Matrix ans;
    for(int i = 0; i < n; ++i)
        for(int j = 0; j < n; ++j)
            ans.v[i][j] = (i == j);
    while (p)
    {
        if (p & 1)
            ans = ans * C;
        C = C * C;
        p >>= 1;
    }
    return ans;
}
Matrix operator + (Matrix a, Matrix b)
{
    Matrix c;
    memset(c.v, 0, sizeof(c.v));
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
                c.v[i][j] = (a.v[i][j] + b.v[i][j]) % mod;
    return c;
}
void output( Matrix a)
{
    REP(i,0,n)
        REP(j,0,n)
            printf("%d%c",a.v[i][j],j<n-1?' ':'\n');
}
Matrix E,A;
Matrix dfs(int r)
{
    if(r==1) return A;
    if( r%2==0 )
    {
        return (matrix_pow(A,r/2) + E ) * dfs(r/2);
    }
    else
    {
        return (matrix_pow(A,r/2) + E) * dfs(r/2) + matrix_pow(A,r);
    }
}

int main()
{
    while(~scanf("%d%d%d",&n,&k,&mod))
    {
        Matrix ans;
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < n; ++j)
                E.v[i][j] = (i==j);
        REP(i,0,n) REP(j,0,n) scan(A.v[i][j]);
        ans = dfs(k);
        output(ans);
    }
}