不久前做过POJ3070,所以知道这题要用矩阵快速幂优化,但是这个题的递推公式中有一项?p/n?,场上就不会了。。。
下来才知道要用分块矩阵快速幂,因为?p/n?最多有2√p块,可以对每一块使用快速幂,复杂度(应该)为lgn*√p。
每一块的范围可以在O(1)的时间内求出,范围为x到min(n,p/(p/x)),具体证明lyd的进阶指南上有。。。
附上代码:
#include<cstdio> #include<algorithm> #include<cstring> using namespace std; typedef long long LL; const int mod=1e9+7; LL a[3][3],f[3]; LL p,n; void mul(LL f[3],LL a[3][3]){ LL c[3]; memset(c,0,sizeof(c)); for(int j=0;j<3;j++) for(int k=0;k<3;k++){ c[j]=(c[j]+f[k]*a[k][j])%mod; } memcpy(f,c,sizeof(c)); } void mulself(LL a[3][3]){ LL c[3][3]; memset(c,0,sizeof(c)); for(int i=0;i<3;i++) for(int j=0;j<3;j++) for(int k=0;k<3;k++) c[i][j]=(c[i][j]+a[i][k]*a[k][j])%mod; memcpy(a,c,sizeof(c)); } void quick_power(int n,LL a[3][3]){ LL b[3][3]; memset(b,0,sizeof(b)); memcpy(b,a,sizeof(a)); for(;n;n>>=1){ if(n&1)mul(f,b); mulself(b); } } int main(){ int T; scanf("%d",&T); while(T--){ memset(a,0,sizeof(a)); a[0][1]=1; f[2]=1; a[2][2]=1; scanf("%lld%lld%lld%lld%lld%lld",&f[1],&f[0],&a[1][0],&a[0][0],&p,&n); if(n==1)printf("%lld\n",f[1]); else if(n==2)printf("%lld\n",f[0]); else{ for(int x=3,gx;x<=n;x=gx+1){ gx=p/x?min(p/(p/x),n):n; a[2][0]=p/x; quick_power(gx-x+1,a); } } printf("%lld\n",f[0]); } }
原文地址:https://www.cnblogs.com/pkgunboat/p/9473573.html
时间: 2024-11-02 22:47:22