以下内容为转载
用二分方法求等比数列前n项和：
原理：
(1)若n==0
(2)若n%2==0
(3)若n%2==1
那么类比等比数列，对于一个矩阵A，求A+A1+A2+.....+A^k
采用分治法，若 k是偶数，原式=（A+A1+.....+A(k/2)）+A^(k/2)*( A+A1+.....+A(k/2) ) = (A^(k/2)+1)*S(k/2).
若k是奇数，原式=（A+A1+.....+A(k/2)）+A^(k/2)*( A+A1+.....+A(k/2) ) + a^k = (A^(k/2)+1)*S(k/2) + A^k.
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
const int MAXN=32;
int mod;
int trie_s;
struct Matrix
{
int arr[MAXN][MAXN];
void init()
{
memset(arr,0,sizeof(arr));
for(int i=1;i<=trie_s;i++)
{
arr[i][i]=1;
}
}
}A,R;
Matrix add(Matrix a,Matrix b)
{
Matrix c;
for(int i=1;i<=trie_s;i++)
{
for(int j=1;j<=trie_s;j++)
{
c.arr[i][j]=a.arr[i][j]+b.arr[i][j];
if(c.arr[i][j]>=mod) c.arr[i][j]%=mod;
}
}
return c;
}
Matrix multi(Matrix a,Matrix b)
{
Matrix c;
memset(c.arr,0,sizeof(c.arr));
for(int i=1;i<=trie_s;i++)
{
for(int k=1;k<=trie_s;k++)
{
if(a.arr[i][k]==0) continue;
for(int j=1;j<=trie_s;j++)
{
c.arr[i][j]+=a.arr[i][k]*b.arr[k][j];
if(c.arr[i][j]>=mod) c.arr[i][j]%=mod;
}
}
}
return c;
}
Matrix pow(Matrix a,int b)
{
Matrix res;
res.init();
while(b)
{
if(b&1) res=multi(res,a);
a=multi(a,a);
b>>=1;
}
return res;
}
Matrix sum(Matrix a,int n)
{
if(n==1) return a;
Matrix tmp;
tmp.init();
tmp=add(tmp,pow(a,n>>1));
tmp=multi(tmp,sum(a,n>>1));
if(n&1) tmp=add(tmp,pow(a,n));
return tmp;
}
void output()
{
for(int i=1;i<=trie_s;i++)
{
for(int j=1;j<=trie_s;j++)
{
printf("%d ",A.arr[i][j]);
}
printf("\n");
}
}
int main()
{
int n;
scanf("%d%d%d",&trie_s,&n,&mod);
for(int i=1;i<=trie_s;i++)
{
for(int j=1;j<=trie_s;j++)
scanf("%d",&A.arr[i][j]);
}
A=sum(A,n);
output();
return 0;
}
还有一种更快的方法，便是构造矩阵。可惜不是我自己想出来的。膜拜，矩阵太强大了。
我们要求的矩阵设为A，先构造这样的矩阵
B = | A A |，
| 0 1 |
B^2 = |A^2 A+A^2|
|0 1 |
....
....
B^k = |A^k A+A2+...+Ak| ，
|0 1 |
因此我们只需求出B ^k，然后取其右上角的n*n的矩阵便是答案
#include<cstdio>
#include<cstring>
using namespace std;
const int MAXN=62;
int mod;
int trie_s;
struct Matrix
{
int arr[MAXN][MAXN];
void init()
{
memset(arr,0,sizeof(arr));
for(int i=1;i<=trie_s;i++)
{
arr[i][i]=1;
}
}
};
Matrix add(Matrix a,Matrix b)
{
Matrix c;
for(int i=1;i<=trie_s;i++)
{
for(int j=1;j<=trie_s;j++)
{
c.arr[i][j]=a.arr[i][j]+b.arr[i][j];
if(c.arr[i][j]>=mod) c.arr[i][j]%=mod;
}
}
return c;
}
Matrix multi(Matrix a,Matrix b)
{
Matrix c;
memset(c.arr,0,sizeof(c.arr));
for(int i=1;i<=trie_s;i++)
{
for(int k=1;k<=trie_s;k++)
{
if(a.arr[i][k]==0) continue;
for(int j=1;j<=trie_s;j++)
{
c.arr[i][j]+=a.arr[i][k]*b.arr[k][j];
if(c.arr[i][j]>=mod) c.arr[i][j]%=mod;
}
}
}
return c;
}
Matrix pow(Matrix a,int b)
{
Matrix res;
res.init();
while(b)
{
if(b&1) res=multi(res,a);
a=multi(a,a);
b>>=1;
}
return res;
}
Matrix sum(Matrix a,int n)
{
if(n==1) return a;
Matrix tmp;
tmp.init();
tmp=add(tmp,pow(a,n>>1));
tmp=multi(tmp,sum(a,n>>1));
if(n&1) tmp=add(tmp,pow(a,n));
return tmp;
}
void output(Matrix A)
{
for(int i=1;i<=trie_s;i++)
{
for(int j=1,j1=trie_s+1;j<=trie_s;j++,j1++)
{
printf("%d ",A.arr[i][j1]);
}
printf("\n");
}
}
int main()
{
int n;
scanf("%d%d%d",&trie_s,&n,&mod);
Matrix A;
for(int i=1;i<=trie_s;i++)
{
for(int j=1;j<=trie_s;j++)
scanf("%d",&A.arr[i][j]);
}
Matrix res;
for(int i=1;i<=trie_s;i++)
{
for(int j=1;j<=trie_s;j++)
{
res.arr[i][j]=A.arr[i][j];
}
}
for(int i=1;i<=trie_s;i++)
{
for(int j=1,j1=trie_s+1;j<=trie_s;j++,j1++)
{
res.arr[i][j1]=A.arr[i][j];
}
}
for(int i1=trie_s+1,i=1;i<=trie_s;i++,i1++)
{
for(int j=1;j<=trie_s;j++)
{
res.arr[i1][j]=0;
}
}
for(int i1=trie_s+1,i=1;i<=trie_s;i++,i1++)
{
for(int j1=trie_s+1,j=1;j<=trie_s;j++,j1++)
{
if(i1!=j1) res.arr[i1][j1]=0;
else res.arr[i1][j1]=1;
}
}
trie_s*=2;
res=pow(res,n);
trie_s/=2;
output(res);
return 0;
}
