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

⎧⎩⎨⎪⎪⎪⎪⎪⎪F1F2Fn===ABC⋅Fn−2+D⋅Fn−1+⌊Pn⌋

  Your job is simple, for each task, you should output Fn module 109+7.

 

Input

The first line has only one integer T, indicates the number of tasks.

Then, for the next T lines, each line consists of 6 integers, A , B, C, D, P, n.

1≤T≤200≤A,B,C,D≤1091≤P,n≤109

 

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));

　　

#include <bits/stdc++.h>
using namespace std;
#define LL long long 
LL mod=1e9+7;
LL f[41000],T,A,B,C,D,P,N;
struct matrix{
    LL a[3][3];
    matrix(){
        memset(a,0,sizeof(a));
    }
    matrix operator*(matrix &tmp){
        matrix ans;
        for(int i=0;i<3;++i){
            for(int j=0;j<3;++j){
                for(int k=0;k<3;++k){
                    ans.a[i][j]=(ans.a[i][j]+a[i][k]*tmp.a[k][j]);
                //    ans.a[i][j]%=mod;
                }
                ans.a[i][j]%=mod;
            }
        }
        return ans;
    }
    
    void show(){
        puts("---");
        for(int i=0;i<3;++i){
            for(int j=0;j<3;++j){
                cout<<a[i][j]<<' ';
            }cout<<endl;
        }
        puts("---");
    }
}X,U;
matrix qpow(matrix A,int b){
    matrix ans=U;
    while(b){
        if(b&1)ans=ans*A;
        A=A*A;
        b>>=1;
    }
    return ans;
}
LL solve(LL i){
    LL key=P/i;
    LL l=i,r=N;
    while(l<r){
        int mid=r-(r-l)/2;
        //cout<<mid<<' '<<P/mid<<endl;
        if(P/mid==key) l=mid;
        else if(P/mid<key){
            r=mid-1;
        }
        else{
            l=mid+1;
        }
    }
    return l;
}
int main(){
    U.a[0][0]=U.a[1][1]=U.a[2][2]=1;
    scanf("%lld",&T);
    while(T--){
        scanf("%lld%lld%lld%lld%lld%lld",&A,&B,&C,&D,&P,&N);
        
        f[1]=A%mod;
        f[2]=B%mod;
        if(N<=2){
            cout<<f[N]<<endl;
            continue;
        }
        //for(int i=3;i<=N;++i) f[i]=(C*f[i-2]%mod+D*f[i-1]%mod+(P/i))%mod;
        memset(X.a,0,sizeof(X.a));
        X.a[0][0]=D,X.a[0][1]=1;
        X.a[1][0]=C,X.a[2][2]=1;
        LL f1=f[2],f2=f[1];
        for(LL i=3;i<=N;){
            LL j=solve(i);
            //cout<<i<<' '<<j<<endl;
            X.a[2][0]=P/i;
            matrix ans=qpow(X,j-i+1);
            LL _f1=(f1*ans.a[0][0]%mod+f2*ans.a[1][0]%mod+ans.a[2][0])%mod;
            LL _f2=(f1*ans.a[0][1]+f2*ans.a[1][1]+ans.a[2][1])%mod;
            
            f1=_f1;
            f2=_f2;
            //cout<<j<<' '<<f1<<' '<<f2<<' '<<f[j]<<' '<<f[j-1]<<endl;
            i=j+1;
        }
        cout<<f1<<endl;
        //printf("%lld %lld
",f1,f[N]);
        
    }    
    return 0;
}
/*
2
3 3 2 1 3 5
3 2 2 2 1 4

*/