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

题目连接

题意：求 A + A2 + A3 + … + Ak. 

数据范围：（n (n ≤ 30), k (k ≤ 109) and m (mod< 104)

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

假设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

;
};
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);
}
}