HDU-6395-矩阵快速幂

Sequence

Time Limit: 4000/2000 MS (Java/Others) Memory Limit: 262144/262144 K (Java/Others)
Total Submission(s): 52 Accepted Submission(s): 13

Problem Description

Let us define a sequence as below

⎧⎩⎨⎪⎪⎪⎪⎪⎪ F 1 F 2 F n = = = A B C \cdot F n - 2 +

Input

The first line has only one integer

Sample Input

2 3 3 2 1 3 5 3 2 2 2 1 4

Sample Output

36 24

Source

2018 Multi-University Training Contest 7

　　对n进行分段，每一段内的 floor(P/ni) 都是一定的，这样就可以对每一段直接跑矩阵快速幂了。分段的时候应该可以O(1)的，比赛时候

没想到直接上的二分= =。

　　(把第78行求边界换成这句代码就是O(1)了: LL j=P/i==0?N:min(N,P/(P/i));

  1 #include <bits/stdc++.h>
  2 using namespace std;
  3 #define LL long long 
  4 LL mod=1e9+7;
  5 LL f[41000],T,A,B,C,D,P,N;
  6 struct matrix{
  7     LL a[3][3];
  8     matrix(){
  9         memset(a,0,sizeof(a));
 10     }
 11     matrix operator*(matrix &tmp){
 12         matrix ans;
 13         for(int i=0;i<3;++i){
 14             for(int j=0;j<3;++j){
 15                 for(int k=0;k<3;++k){
 16                     ans.a[i][j]=(ans.a[i][j]+a[i][k]*tmp.a[k][j]);
 17                 //    ans.a[i][j]%=mod;
 18                 }
 19                 ans.a[i][j]%=mod;
 20             }
 21         }
 22         return ans;
 23     }
 24     
 25     void show(){
 26         puts("---");
 27         for(int i=0;i<3;++i){
 28             for(int j=0;j<3;++j){
 29                 cout<<a[i][j]<<' ';
 30             }cout<<endl;
 31         }
 32         puts("---");
 33     }
 34 }X,U;
 35 matrix qpow(matrix A,int b){
 36     matrix ans=U;
 37     while(b){
 38         if(b&1)ans=ans*A;
 39         A=A*A;
 40         b>>=1;
 41     }
 42     return ans;
 43 }
 44 LL solve(LL i){
 45     LL key=P/i;
 46     LL l=i,r=N;
 47     while(l<r){
 48         int mid=r-(r-l)/2;
 49         //cout<<mid<<' '<<P/mid<<endl;
 50         if(P/mid==key) l=mid;
 51         else if(P/mid<key){
 52             r=mid-1;
 53         }
 54         else{
 55             l=mid+1;
 56         }
 57     }
 58     return l;
 59 }
 60 int main(){
 61     U.a[0][0]=U.a[1][1]=U.a[2][2]=1;
 62     scanf("%lld",&T);
 63     while(T--){
 64         scanf("%lld%lld%lld%lld%lld%lld",&A,&B,&C,&D,&P,&N);
 65         
 66         f[1]=A%mod;
 67         f[2]=B%mod;
 68         if(N<=2){
 69             cout<<f[N]<<endl;
 70             continue;
 71         }
 72         //for(int i=3;i<=N;++i) f[i]=(C*f[i-2]%mod+D*f[i-1]%mod+(P/i))%mod;
 73         memset(X.a,0,sizeof(X.a));
 74         X.a[0][0]=D,X.a[0][1]=1;
 75         X.a[1][0]=C,X.a[2][2]=1;
 76         LL f1=f[2],f2=f[1];
 77         for(LL i=3;i<=N;){
 78             LL j=solve(i);
 79             //cout<<i<<' '<<j<<endl;
 80             X.a[2][0]=P/i;
 81             matrix ans=qpow(X,j-i+1);
 82             LL _f1=(f1*ans.a[0][0]%mod+f2*ans.a[1][0]%mod+ans.a[2][0])%mod;
 83             LL _f2=(f1*ans.a[0][1]+f2*ans.a[1][1]+ans.a[2][1])%mod;
 84             
 85             f1=_f1;
 86             f2=_f2;
 87             //cout<<j<<' '<<f1<<' '<<f2<<' '<<f[j]<<' '<<f[j-1]<<endl;
 88             i=j+1;
 89         }
 90         cout<<f1<<endl;
 91         //printf("%lld %lld
",f1,f[N]);
 92         
 93     }    
 94     return 0;
 95 }
 96 /*
 97 2
 98 3 3 2 1 3 5
 99 3 2 2 2 1 4
100 
101 */