hdu 6395 整数分块+矩阵快速幂

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
typedef struct node
{
    ll mat[4][4];
}matrix;
matrix mul(matrix a,matrix b)
{
    matrix res;
    for(int i=0;i<4;i++)
        for(int j=0;j<4;j++)
            res.mat[i][j]=0;
    for(int i=0;i<4;i++)
        for(int j=0;j<4;j++)
            for(int k=0;k<4;k++)
                res.mat[i][j]=(res.mat[i][j]+a.mat[i][k]*b.mat[k][j]%mod)%mod;
    return res;
}
matrix quick_mod(matrix n,ll m)
{
    matrix res;
    for(int i=0;i<4;i++)
    {
        for(int j=0;j<4;j++)
            res.mat[i][j]=0;
        res.mat[i][i]=1;
    }
    while(m)
    {
        if(m&1)
            res=mul(res,n);
        n=mul(n,n);
        m>>=1;
    }
    return res;
}
int main()
{
    ll t,a,b,c,d,p,n,last,tpp,ans;
    scanf("%lld",&t);
    while(t--&&scanf("%lld%lld%lld%lld%lld%lld",&a,&b,&c,&d,&p,&n)!=EOF)
    {
        matrix res;
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4;j++)
                res.mat[i][j]=0;
            res.mat[i][i]=1;
        }
        matrix tp;
        tp.mat[0][0]=d;
        tp.mat[1][0]=c;
        tp.mat[0][1]=tp.mat[3][3]=tp.mat[2][0]=tp.mat[2][2]=1;
        tp.mat[0][2]=tp.mat[1][1]=tp.mat[1][2]=tp.mat[1][3]=tp.mat[2][1]=tp.mat[2][3]=tp.mat[3][1]=tp.mat[0][3]=0;
        if(n<=2)
        {
            if(n==1)
                printf("%d\n",a);
            else if(n==2)
                printf("%d\n",a);
            continue;
        }
        tpp=n;
        if(n>p)
        {
            if(p>=3)
            {
                tp.mat[3][0]=tp.mat[3][2]=-1;
                res=mul(res,tp);
                tp.mat[3][0]=tp.mat[3][2]=0;
                res=mul(res,quick_mod(tp,n-p-1));
                n=p;
            }
            else
            {
                tp.mat[3][0]=tp.mat[3][2]=0;
                res=mul(res,quick_mod(tp,n-2));
                n=2;
            }
        }
        if(p<5)
        {
            ll fn,fn1=b,fn2=a;
            for(int i=3;i<=n&&i<=tpp;i++)
            {
                fn=(c*fn2%mod+d*fn1%mod+p/i)%mod;
                fn2=fn1;
                fn1=fn;
            }
            if(tpp>p)
            {
                ans=(fn1*res.mat[0][0]%mod+fn2*res.mat[1][0]%mod+p/n*res.mat[2][0]+res.mat[3][0])%mod;
                printf("%lld\n",ans);
            }
            else
                printf("%lld\n",fn);
            continue;
        }
        matrix ppt=res;
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4;j++)
                res.mat[i][j]=0;
            res.mat[i][i]=1;
        }
        last=p/2;
        tpp=(int)(sqrt((p-1)*1.0)+1);
        for(int j=3;j<=tpp&&j<=n;j++)
        {
            tp.mat[3][0]=tp.mat[3][2]=p/j-last;
            last=p/j;
            res=mul(res,tp);
        }
        if(n>tpp)
        {
            ll now=tpp;
            now--;
            while(last)
            {
                tpp=p/last;
                tp.mat[3][0]=tp.mat[3][2]=0;
                if(tpp>=n)
                {
                    res=mul(res,quick_mod(tp,n-now-1));
                    break;
                }
                else
                    res=mul(res,quick_mod(tp,tpp-now-1));
                now=tpp;
                tp.mat[3][0]=tp.mat[3][2]=-1;
                res=mul(res,tp);
                last--;
            }
        }
        res=mul(res,ppt);
        ans=((b*res.mat[0][0]%mod+a*res.mat[1][0]%mod+p/2*res.mat[2][0]%mod+res.mat[3][0])%mod+mod)%mod;
        printf("%lld\n",ans);
    }
    return 0;
}