Sequence(整除分块+矩阵快速幂)

2018-08-14

Sequence

题意

定义：
$$
F_1=A
$$

$$
F_2=B
$$

$$
F_n=D\cdot{}F_{n-1}+C\cdot{}F_{n-2}+\lfloor\frac{P}{n}\rfloor
$$

求Fn，Mod1e9+7。

思路

一看到这个向下取整！

多么亲切啊！

整除分块啊！

……

一开始怎么都想不出来怎么把分块用进去，队友一波机智，先分块，然后矩阵快速幂，哇，秒啊，然后两个小时过去了……

由于莫比乌斯中我们经常会用到一个分块的概念：
$$
\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor
$$

for(int l=1,r;l<=n;l=r+1)
{
    r=n/(n/l);
    ans+=(r-l+1)*(n/l);
}

那么我们就想，因为矩阵快速幂要求的是尾项是个常数，而对于这个整除来说，在某些区间内它的值是一定的，因此我们就可以先把他们划分为好几块，然后对于每个区间进行矩阵快速幂即可，有些细节注意，比如说次数不足3次时，要手写一波判断，因为没有需要的Fn-1，Fn-2，比如n小于p的时候，需要在循环内直接结束，比如n大于p的时候，分块结束后再进行裸的矩阵快速幂。

看！代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD = 1e9+7;
struct mat
{
    long long a[3][3];
}m,r;
mat mult(mat x,mat y)
{
    mat res={0};
    int i,j,k;
    for(i=0;i<3;i++)
    for(j=0;j<3;j++)
    for(k=0;k<3;k++)
        res.a[i][j]=(res.a[i][j]+(x.a[i][k]*y.a[k][j])%MOD)%MOD;
    return res;
}
mat PowerMod(mat x,long long n)
{
    mat ans={0};
    int i,j;
    for(i=0;i<3;i++)
    ans.a[i][i]=1;
    while(n>0)
    {
        if(n%2==1)
        ans=mult(ans,x);
        x=mult(x,x);
        n=n/2;
    }
    return ans;
}
int main()
{
    int t;
    ll a,b,c,d,p,n;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld%lld%lld%lld",&a,&b,&c,&d,&p,&n);
        if(n==1) { printf("%lld\n",a);continue;}
        else if(n==2) { printf("%lld\n",b);continue;}
        m.a[0][2]=m.a[1][1]=m.a[2][1]=m.a[2][0]=0;
        m.a[0][0]=d;m.a[0][1]=c;
        m.a[1][0]=m.a[2][2]=1;
        ll pre=b,ppre=a;
        for(int L=3,R;L<=p;L=R+1)
        {

            R=p/(p/L); 
            int zs=min(1ll*R-L+1,n-L+1);
            if(zs<=3)
            {
                for(int i=1;i<=zs;i++)
                {
                    ll tmp=pre;
                    pre=((pre*d+ppre*c)%MOD+p/L)%MOD;
                    ppre=tmp;
                }
                if(n>=L&&n<=R) break;
                continue;
            }
            m.a[0][2]=p/L;
            r=PowerMod(m,zs);
            ll p1=(r.a[0][0]*pre+r.a[0][1]*ppre+r.a[0][2])%MOD;
            ll p2=(r.a[1][0]*pre+r.a[1][1]*ppre+r.a[1][2])%MOD;
            pre=p1;ppre=p2;
            if(n>=L&&n<=R) break;
        }
        if(n>p)
        {
            m.a[0][2]=0;
            r=PowerMod(m,n-max(2ll,p));
            ll p1=(r.a[0][0]*pre+r.a[0][1]*ppre+r.a[0][2])%MOD;
            pre=p1;
        }
        printf("%lld\n",pre);
    }
    return 0;
}