Sequence

Problem Description
Let us define a sequence as below

f1=A
f2=B
fn=C*fn-2+D*fn-1+[p/n]

  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 , BCDPn.

1T200A,B,C,D1091P,n109
 
Sample Input
2
3 3 2 1 3 5
3 2 2 2 1 4
 
Sample Output
36 24
 
题意:题目给出ABCDPn,第一项是A,第二项是B,然后还有个递推式,问第n项是多少
 
思路:如果我们按照他的递推式去推得答案的话,n的范围是1e9肯定会超时,但是我们又必须要用到这个式子,我们就想有没有加快的方法,其实做多了题会发现这是矩阵快速幂的形式
矩阵快速幂就是把你原有的递推式再加快执行,
但是 fn=C*fn-2+D*fn-1+[p/n]
其中[p/n]是个会变化的值,我们的矩阵里面不好处理
如果是一个常数C1的话
我们的关系矩阵就可以写为 
fn  D  C  1    fn-1
fn-1=  1   0  0 *    fn-2
C1      0     0  1     C1 
 利用了矩阵乘法
 
然后我们再来讨论[p/n]
[p/n]因为是整除,那么他的一段内的值肯定是相同的
这个我们可以利用数论中的除法分块得出 [p/n]   n属于1-p   他的整除值也只会有 O(sqrt(p))种值
我们对每一块使用矩阵快速幂,那样这个整除数又是常数就可以利用上面这个关系矩阵,由于只有O(sqrt(p))种值,复杂度上也不会有事
所以我们总的解法就是 除法分块+矩阵快速幂
 
除法分块:i-p/(p/i) 这一段是一个块
因为他们的值都是 p/i  ,最小的是i,因为是余数最多的情况,我们利用最基本的除法原理   a/b=c   ->    a/c=b   ,我们把整除值当作除数 ,除出来的自然是下标,又因为是整除,所以
算出来的下标也是余数最小的情况   
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<ctime>
#include<algorithm>
using namespace std;
const long long mod=1e9+7;
struct jz//结构体写法的矩阵快速幂
{
    long long num[3][3];
    jz() { memset(num,0,sizeof(num)); }
    jz operator*(const jz &P)const
    {
        jz ans;
        for(int k=0;k<3;k++)
            for(int i=0;i<3;i++)
                for(int j=0;j<3;j++)
                    ans.num[i][j]=(ans.num[i][j]+num[i][k]*P.num[k][j]%mod)%mod;
        return ans;
    }
}COE,ans,unit;
int T_T;
long long A,B,C,D,P,n;
jz pOw(jz X,long long m)//矩阵快速幂
{
    jz ans;
    for(ans=unit;m;m>>=1,X=X*X)
        if(m&1)
            ans=ans*X;
    return ans;
}
void init(long long A,long long B,long long C,long long D,long long x)//更新关系矩阵
{
    COE.num[0][0]=0;
    COE.num[0][1]=1;
    COE.num[0][2]=0;
    COE.num[1][0]=C;
    COE.num[1][1]=D;
    COE.num[1][2]=x; // this element need to be changed each step.
    COE.num[2][0]=0;
    COE.num[2][1]=0;
    COE.num[2][2]=1;
    return;
}
int main()
{
    for(int i=0;i<3;i++) unit.num[i][i]=1;
    scanf("%d",&T_T);
    while(T_T--)
    {
        scanf("%lld%lld%lld%lld%lld%lld",&A,&B,&C,&D,&P,&n);
        if(n==1) printf("%lld\n",A);
        else if(n<P)
        {
            ans.num[0][0]=A;
            ans.num[1][0]=B;
            ans.num[2][0]=1;
            for(long long i=3;i<=n;i=P/(P/i)+1)//除法分块
            {
                init(A,B,C,D,P/i);
                if(n<=P/(P/i)) COE=pOw(COE,n-i+1);
                else COE=pOw(COE,P/(P/i)+1-i);
                ans=COE*ans;
            }
            printf("%lld\n",ans.num[1][0]);
        }
        else if(P<=n)
        {
            ans.num[0][0]=A;
            ans.num[1][0]=B;
            ans.num[2][0]=1;
            for(long long i=3;i<=P;i=P/(P/i)+1)//除法分块
            {
                init(A,B,C,D,P/i);
                COE=pOw(COE,P/(P/i)+1-i);
                ans=COE*ans;
            }
            init(A,B,C,D,0);
            COE.num[1][2]=0;
            COE=pOw(COE,n-max(P,2LL));//多余的一段的整除值都是0
            ans=COE*ans;
            printf("%lld\n",ans.num[1][0]);
        }
    }
    fclose(stdin);
    fclose(stdout);
    return 0;
}