hdu6395-分段矩阵快速幂

矩阵快速幂加整除分块

#include <iostream>
#include <cstdio>
#include <cstring>
#include <map>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
struct Matrix
{
    ll wm[3][3];
};
Matrix E;
void Copy(Matrix *a,Matrix *b)
{
    for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
    a->wm[i][j]=b->wm[i][j];
}
void Matrix_mul(Matrix *a,Matrix *b,Matrix *c)
{
    Matrix t;
    memset(t.wm,0,sizeof(t.wm));
    for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
    for(int k=0;k<3;k++) t.wm[i][j]+=a->wm[i][k]*b->wm[k][j],t.wm[i][j]%=mod;
    Copy(c,&t);
}
void Matrix_pow(Matrix *a,ll n,Matrix *c)
{
    Matrix base,ans;
    Copy(&ans,&E),Copy(&base,a);
    while(n)
    {
        if(n&1) Matrix_mul(&ans,&base,&ans);
        Matrix_mul(&base,&base,&base);
        n>>=1;
    }
    Copy(c,&ans);
}

void init(Matrix *Base,ll c,ll d)
{
    for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
    if(i==j) E.wm[i][j]=1;
    else E.wm[i][j]=0;
    Base->wm[0][0]=d,Base->wm[0][1]=c,Base->wm[0][2]=1;
    Base->wm[1][0]=1,Base->wm[1][1]=0,Base->wm[1][2]=0;
    Base->wm[2][0]=0,Base->wm[2][1]=0,Base->wm[2][2]=1;
}
void compute_F(Matrix *Base,ll n,Matrix *pre)
{
    if(n==0) {
            memset(pre->wm,0,sizeof(pre->wm));
    return;
    }
    if(n==1) {
    Copy(pre,Base);
    return ;
    }
    Matrix temp;
    Matrix_pow(Base,n-1,&temp);
    Copy(pre,&temp);
}
inline ll max(ll a,ll b)
{
    if(a>b) return a;
    return b;
}
inline ll min(ll a,ll b)
{
    if(a>b) return b;
    return a;
}
int main()
{
    int t;
    ll a,b,c,d,p,n,now1,now2,l,r;
    Matrix Base,temp;
  //  freopen("in.txt","r",stdin);
  //  freopen("out.txt","w",stdout);
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld%lld%lld%lld",&a,&b,&c,&d,&p,&n);
        if(n==1)
            b=a;
        else if(n==2)
            b=b;
        else
        {
            init(&Base,c,d);
            Matrix pre;
            for(l=3;l<=min(p,n);l=r+1)
            {
                r=min(n,p/(p/l));
                Matrix_pow(&Base,r-l+1,&pre);
                now1=(pre.wm[0][0]*b%mod+pre.wm[0][1]*a%mod+pre.wm[0][2]*(p/l)%mod)%mod;
                now2=(pre.wm[1][0]*b%mod+pre.wm[1][1]*a%mod+pre.wm[1][2]*(p/l)%mod)%mod;
                b=now1,a=now2;
            }
            p=max((ll)2,p);
            if(n>p)
            {
                Matrix_pow(&Base,n-p,&pre);
                now1=(pre.wm[0][0]*b%mod+pre.wm[0][1]*a%mod)%mod;
                b=now1;
            }
        }
        printf("%lld\n",b);
    }
    return 0;
}