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

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

struct jz//结构体写法的矩阵快速幂

{

    long long num[][];

    jz() { memset(num,,sizeof(num)); }

    jz operator*(const jz &P)const

    {

        jz ans;

        for(int k=;k<;k++)

            for(int i=;i<;i++)

                for(int j=;j<;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>>=,X=X*X)

        if(m&)

            ans=ans*X;

    return ans;

}

void init(long long A,long long B,long long C,long long D,long long x)//更新关系矩阵

{

    COE.num[][]=;

    COE.num[][]=;

    COE.num[][]=;

    COE.num[][]=C;

    COE.num[][]=D;

    COE.num[][]=x; // this element need to be changed each step.

    COE.num[][]=;

    COE.num[][]=;

    COE.num[][]=;

    return;

}

int main()

{

    for(int i=;i<;i++) unit.num[i][i]=;

    scanf("%d",&T_T);

    while(T_T--)

    {

        scanf("%lld%lld%lld%lld%lld%lld",&A,&B,&C,&D,&P,&n);

        if(n==) printf("%lld\n",A);

        else if(n<P)

        {

            ans.num[][]=A;

            ans.num[][]=B;

            ans.num[][]=;

            for(long long i=;i<=n;i=P/(P/i)+)//除法分块

            {

                init(A,B,C,D,P/i);

                if(n<=P/(P/i)) COE=pOw(COE,n-i+);

                else COE=pOw(COE,P/(P/i)+-i);

                ans=COE*ans;

            }

            printf("%lld\n",ans.num[][]);

        }

        else if(P<=n)

        {

            ans.num[][]=A;

            ans.num[][]=B;

            ans.num[][]=;

            for(long long i=;i<=P;i=P/(P/i)+)//除法分块

            {

                init(A,B,C,D,P/i);

                COE=pOw(COE,P/(P/i)+-i);

                ans=COE*ans;

            }

            init(A,B,C,D,);

            COE.num[][]=;

            COE=pOw(COE,n-max(P,2LL));//多余的一段的整除值都是0

            ans=COE*ans;

            printf("%lld\n",ans.num[][]);

        }

    }

    fclose(stdin);

    fclose(stdout);

    return ;

}

秒客网

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

Sequence