[TopCoder11557]MatrixPower

vjudge

description

你有一个\(n \times n\)的矩阵\(A\)，下标从\(0\)开始，其中\(A_{i,j}=di + q^j\)。
给你\(d,q,n,k,s,t\)，求\((A^k)_{s,t}\)的值模\(10^9+7\)。
\(n\le10^4,d,q,k\le10^9,0\le s,t <n\)

sol

先写出一个极其丑陋的柿子：

\[(A^k)_{s,t}=\sum_{x_1}\sum_{x_2}...\sum_{x_{k-1}}(ds+q^{x_1})(dx_1+q^{x_2})...(dx_{k-1}+q^t) \]

等于你从\(s\)出发走\(k\)次走到\(t\)，那么中间就有\(k-1\)个位置可以任选，所以才会有前面的\(k-1\)个\(\sum\)。
把后面的括号暴力展开，一共有\(2^k\)项。然而还是做不了。
考虑每个\(x_i\)。它出现在答案中的形式只有\(4\)种：\(1,dx_i,q_{x_i},dx_iq^{x_i}\)。这取决于这个\(x_i\)出现的两个相邻的括号内取的是哪一个。
预处理出\(sd=\sum_{i=0}^{n-1}di,sq=\sum_{i=0}^{n-1}q^i,ss=\sum_{i=0}^{n-1}diq^i\)。
设\(f_{i,0/1}\)表示考虑了前\(i\)个括号，第\(i\)个括号取的是前一项/后一项时，只计算前\(i-1\)个\(x_i\)时的答案。显然有转移式：

\[f_{i+1,0}=f_{i,0}\times sd+f_{i,1}\times ss\\f_{i+1,1}=f_{i,0}\times n+f_{i,1}\times sq \]

初值\(f_{1,0}=ds,f_{1,1}=1\)。最后一次转移特判一下（因为最后一个\(x\)必须强制为\(t\)）。
这样就可以做\(O(k)\)的\(dp\)了。考虑到每次的转移都是一样的，所以直接矩乘就好了。

code

为什么没有人告诉我，定义在\(\mbox{class}\)里面的变量的初值不是\(0\)。
害我浪费了宝贵的\(\mbox{20min}\)才发现这个错误\(\mbox{qwq}\)。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
using namespace std;
int gi(){
	int x=0,w=1;char ch=getchar();
	while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
	if (ch=='-') w=0,ch=getchar();
	while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
	return w?x:-x;
}
const int N = 10005;
const int mod = 1e9+7;
int fastpow(int a,int b){
	int s=1;
	while (b) {if (b&1) s=1ll*s*a%mod;a=1ll*a*a%mod;b>>=1;}
	return s;
}
struct matrix{
	int a[2][2];
	matrix(){memset(a,0,sizeof(a));}
	matrix(int q,int w,int d,int b){
		a[0][0]=q;a[0][1]=w;a[1][0]=d;a[1][1]=b;
	}
	int* operator [](int x){return a[x];}
	matrix operator * (matrix &b){
		matrix c;
		for (int i=0;i<2;++i)
			for (int j=0;j<2;++j)
				for (int k=0;k<2;++k)
					c[i][k]=(1ll*a[i][j]*b[j][k]+c[i][k])%mod;
		return c;
	}
}S,T,A,B;
class MatrixPower{
public:
	int m,sd,sq,ss;vector<int>ans;
	vector<int> getElements(int d,int q,int n,int k,vector<int>r,vector<int>c){
		m=r.size();ans.resize(m);
		if (k==1){
			for (int i=0;i<m;++i) ans[i]=(1ll*d*r[i]+fastpow(q,c[i]))%mod;
			return ans;
		}
		k-=2;sd=sq=ss=0;
		for (int i=0;i<n;++i) sd=(sd+1ll*d*i)%mod;
		for (int i=0;i<n;++i) sq=(sq+fastpow(q,i))%mod;
		for (int i=0;i<n;++i) ss=(ss+1ll*d*i%mod*fastpow(q,i))%mod;
		T=matrix(sd,n,ss,sq);S=matrix(1,0,0,1);
		while (k) {if (k&1) S=S*T;T=T*T;k>>=1;}
		for (int i=0;i<m;++i){
			A=matrix(1ll*d*r[i]%mod,1,0,0);
			B=matrix(sd,1ll*n*fastpow(q,c[i])%mod,ss,1ll*sq*fastpow(q,c[i])%mod);
			A=A*S*B;
			ans[i]=(1ll*A[0][0]+A[0][1]+A[1][0]+A[1][1])%mod;
		}
		return ans;
	}
};