杭电多校第七场 1010 Sequence（除法分块+矩阵快速幂）

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

 

题意：题目给出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
",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
",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
",ans.num[1][0]);
        }
    }
    fclose(stdin);
    fclose(stdout);
    return 0;
}