题面
(T) 组测试数据,每次给定 (n,k),求((F(i)) 为斐波那契数列第 (i) 项):
[sum_{1le x_ile n(1le ile k)}F(sum x_i-k+1) ]
数据范围:(1le Tle 100),(1le n,kle 10^9)。
蒟蒻语
小蒟蒻上午的训练赛头很晕,于是就一直在干这题。
推了 (10) 页草稿纸,推出了很多巧妙的东西(
可惜,天有绝蒻之路,小蒟蒻最后没做出来,幸得爆零。
赛后神仙教了蒟蒻做法,太巧妙了,比出题人的题解强好多倍。
蒟蒻解
把所有 (x_i) 都 (-1),开始推式:
[egin{split}
&sum_{0le x_ile n-1}F(sum x_i+1)\
=&sum_{0le x_ile n-1}
egin{bmatrix}
1&0\
end{bmatrix}
egin{bmatrix}
1&1\
1&0\
end{bmatrix}
^{sum x_i}
egin{bmatrix}
1\
0\
end{bmatrix}\
=&
egin{bmatrix}
1&0\
end{bmatrix}
left(
sum_{0le x_ile n-1}
egin{bmatrix}
1&1\
1&0\
end{bmatrix}
^{sum x_i}
ight)
egin{bmatrix}
1\
0\
end{bmatrix}\
=&
egin{bmatrix}
1&0\
end{bmatrix}
left(
sum_{i=0}^{n-1}
egin{bmatrix}
1&1\
1&0\
end{bmatrix}
^i
ight)^k
egin{bmatrix}
1\
0\
end{bmatrix}
end{split}
]
现在唯一的问题是,如何 (Theta(log n)) 求出 (sumlimits_{i=0}^{n-1}egin{bmatrix}1&1\1&0\end{bmatrix}^i)?
当 (i=0) 的时候,(egin{bmatrix}1&1\1&0\end{bmatrix}^i=E),(E) 是单位矩阵,(E_{i,i}=1),其他位置为 (0)。
神仙的做法是倍增,维护:
[sm_k=sum_{i=0}^{2^k-1}egin{bmatrix}1&1\1&0\end{bmatrix}^i
]
然后这个东西是可以 (Theta(log n)) 递推的,于是 (sumlimits_{i=0}^{n-1}egin{bmatrix}1&1\1&0\end{bmatrix}^i) 就可以 (Theta(log n)) 计算了。
另一种方法是矩阵开两倍空间,做矩阵快速幂求矩阵。
但是无论用哪种方法,求矩阵的 (k) 次幂还是要矩阵快速幂的。
代码
#include <bits/stdc++.h>
using namespace std;
//Start
typedef long long ll;
typedef double db;
#define mp(a,b) make_pair(a,b)
#define x first
#define y second
#define be(a) a.begin()
#define en(a) a.end()
#define sz(a) int((a).size())
#define pb(a) push_back(a)
const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;
//Data
const int mod=1e9+7;
//Matrix
const int M=2;
struct Matrix{
int arr[2][2];
Matrix(){memset(arr,0,sizeof arr);}
Matrix(int x,int y,int z,int w){
arr[0][0]=x,arr[0][1]=y;
arr[1][0]=z,arr[1][1]=w;
}
int*operator[](int x){return arr[x];}
friend Matrix operator+(Matrix a,Matrix b){
Matrix res;
for(int i=0;i<M;i++)
for(int j=0;j<M;j++)
res[i][j]=(a[i][j]+b[i][j])%mod;
return res;
}
friend Matrix operator*(Matrix a,Matrix b){
Matrix res;
for(int k=0;k<M;k++)
for(int i=0;i<M;i++)
for(int j=0;j<M;j++)
(res[i][j]+=(ll)a[i][k]*b[k][j]%mod)%=mod;
return res;
}
};
const Matrix E(1,0,0,1);
Matrix Pow(Matrix a,int x){
Matrix res(E);
for(;x;a=a*a,x>>=1)if(x&1) res=res*a;
return res;
}
void Print(Matrix a){
cout<<"++++++++++++
";
cout<<a[0][0]<<','<<a[0][1]<<'
';
cout<<a[1][0]<<','<<a[1][1]<<'
';
cout<<"------------
";
}
Matrix key(1,0,0,0);
Matrix pa[32],sm[32];
//Main
int main(){
ios::sync_with_stdio(0);
cin.tie(0),cout.tie(0);
sm[0]=E,pa[0]=Matrix(1,1,1,0);
for(int i=1;i<32;i++){
pa[i]=pa[i-1]*pa[i-1];
sm[i]=sm[i-1]+pa[i-1]*sm[i-1];
}
// for(int i=0;i<5;i++){
// cout<<"i="<<i<<'
';
// Print(sm[i]),Print(pa[i]);
// }
int T; cin>>T;
while(T--){
int n,k; cin>>n>>k;
Matrix res;
for(int i=0;i<32;i++)if(n>>i&1)
res=res*pa[i]+sm[i];
res=key*Pow(res,k)*key;
cout<<res[0][0]<<'
';
}
return 0;
}
祝大家学习愉快!