http://poj.org/problem?id=3233
ps转:
用二分方法求等比数列前n项和:即
原理:
(1)若n==0
(2)若n%2==0
(3)若n%2==1
代码如下:
LL sum(LL p,LL n)
{
if(n==0) return 1;
if(n&1) return (1+pow(p,(n>>1)+1))*sum(p,n>>1);
else return (1+pow(p,(n>>1)+1))*sum(p,(n-1)>>1)+pow(p,n>>1);
}
那么类比等比数列,对于一个矩阵A,求A+A^1+A^2+.....+A^k
采用分治法,若 k是偶数,原式=(A+A^1+.....+A^(k/2))+A^(k/2)*( A+A^1+.....+A^(k/2) ) = (A^(k/2)+1)*S(k/2).
若k是奇数,原式=(A+A^1+.....+A^(k/2))+A^(k/2)*( A+A^1+.....+A^(k/2) ) + a^k = (A^(k/2)+1)*S(k/2) + A^k.
我们可以二分等比矩阵中的幂K,每次分成一半求。
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <set>
#include <map>
#include <vector>
#include <math.h>
#include <string.h>
#include <queue>
#include <string>
#include <stdlib.h>
#define LL long long
#define _LL __int64
#define eps 1e-8
#define PI acos(-1.0)
using namespace std;
const int maxn = 35;
int mod,n,k,m;
struct matrix
{
int mat[maxn][maxn];
//初始化为单位矩阵
void init()
{
memset(mat,0,sizeof(mat));
for(int i = 1; i <= maxn; i++)
mat[i][i] = 1;
}
}a,res;
//矩阵相加
matrix matrixAdd(matrix x, matrix y)
{
matrix tmp;
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n; j++)
{
tmp.mat[i][j] = x.mat[i][j] + y.mat[i][j];
if( tmp.mat[i][j] >= mod )
tmp.mat[i][j] %= mod;
}
return tmp;
}
//矩阵相乘
matrix matrixMul(matrix x, matrix y)
{
matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat));
for(int i = 1; i <= n; i++)
{
for(int k = 1; k <= n; k++)
{
if(x.mat[i][k] == 0) continue;
for(int j = 1; j <= n; j++)
{
tmp.mat[i][j] += x.mat[i][k] * y.mat[k][j];
if(tmp.mat[i][j] >= mod)
tmp.mat[i][j] %= mod;
}
}
}
return tmp;
}
//矩阵求幂
matrix Mul(matrix x, int k)
{
matrix tmp;
tmp.init();
while(k)
{
if(k & 1)
tmp = matrixMul(tmp,x);
x = matrixMul(x,x);
k >>= 1;
}
return tmp;
}
// 求A + A^1 + A^2 + ... + A^k
matrix Sum(matrix x, int k)
{
if(k == 1)
return x;
matrix tmp;
tmp.init();
tmp = matrixAdd(tmp,Mul(x,k>>1));
tmp = matrixMul(tmp,Sum(x,k>>1));
if(k&1)
tmp = matrixAdd(tmp,Mul(x,k));
return tmp;
}
void output()
{
for(int i = 1; i <= n; i++)
{
for(int j = 1; j <= n; j++)
{
if(j < n)
printf("%d ",res.mat[i][j]);
else printf("%d\n",res.mat[i][j]);
}
}
}
int main()
{
scanf("%d %d %d",&n,&k,&m);
mod = m;
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n; j++)
scanf("%d",&a.mat[i][j]);
res = Sum(a,k);
output();
return 0;
}
还有一种更快的方法,便是构造矩阵。可惜不是我自己想出来的。膜拜,矩阵太强大了。
我们要求的矩阵设为A,先构造这样的矩阵
B = |A A|,可以计算B^2 = |A^2 A+A^2| .... B^k = |A^k A+A^2+...+A^k| ,
|0 1| |0 1 | |0 1 |
因此我们只需求出B ^k,然后取其右上角的n*n的矩阵便是答案。。
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <set>
#include <map>
#include <vector>
#include <math.h>
#include <string.h>
#include <queue>
#include <string>
#include <stdlib.h>
#define LL long long
#define _LL __int64
#define eps 1e-8
#define PI acos(-1.0)
using namespace std;
const int maxn = 100;
int mod,n,k,m;
struct matrix
{
int mat[maxn][maxn];
//初始化为单位矩阵
void init()
{
memset(mat,0,sizeof(mat));
for(int i = 1; i <= n; i++)
mat[i][i] = 1;
}
}a,res;
//矩阵相加
matrix matrixAdd(matrix x, matrix y)
{
matrix tmp;
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n; j++)
{
tmp.mat[i][j] = x.mat[i][j] + y.mat[i][j];
if( tmp.mat[i][j] >= mod )
tmp.mat[i][j] %= mod;
}
return tmp;
}
//矩阵相乘
matrix matrixMul(matrix x, matrix y)
{
matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat));
for(int i = 1; i <= n; i++)
{
for(int k = 1; k <= n; k++)
{
if(x.mat[i][k] == 0) continue;
for(int j = 1; j <= n; j++)
{
tmp.mat[i][j] += x.mat[i][k] * y.mat[k][j];
if(tmp.mat[i][j] >= mod)
tmp.mat[i][j] %= mod;
}
}
}
return tmp;
}
//矩阵求幂
matrix Mul(matrix x, int k)
{
matrix tmp;
tmp.init();
while(k)
{
if(k & 1)
tmp = matrixMul(tmp,x);
x = matrixMul(x,x);
k >>= 1;
}
return tmp;
}
void output()
{
for(int i = 1; i <= n; i++)
{
for(int j = n+1; j <= n*2; j++)
{
if(j < n*2)
printf("%d ",res.mat[i][j]);
else printf("%d\n",res.mat[i][j]);
}
}
}
int main()
{
while(~scanf("%d %d %d",&n,&k,&m))
{
mod = m;
//扩大A矩阵,看成是4个小矩阵组合而成
for(int i = 1; i <= n; i++)
for(int j = 1; j <= n; j++)
scanf("%d",&a.mat[i][j]);
for(int i = 1; i <= n; i++)
for(int j = n+1; j <= n*2; j++)
a.mat[i][j] = a.mat[i][j-n];
for(int i = n+1; i <= n*2; i++)
for(int j = 1; j <= n; j++)
a.mat[i][j] = 0;
for(int i = n+1; i <= n*2; i++)
for(int j = n+1; j <= n*2; j++)
n = n << 1;
res = Mul(a,k);
n = n >> 1;
output();
}
return 0;
}