HDU-6395-Sequence (矩阵快速密, 加整数分块)

 多校虐我千百遍!

题目很清晰,一看就是矩阵快速幂,但是p/n是个很头疼的问题,比赛的时候想过把所有相同的p/n分一块,但是不会分,赛后看了题解,发现有这么一个神奇的规律:从i到p/(p/i)区间里的数p/i的结果是一样的。这个结论好像在紫书上?

思路:题目给的递推式是一个非线性的关系,所以一开始就想转化成线性关系然后再将转移矩阵求出来,但是这道题没法转化成线性关系直接得到转移矩阵,因此只能根据p/n的值来分段求,因为p/n对于不同的连续的n在一段范围内是一样,因此就可以转化成Fn = C*Fn-2+D*Fn-1+K(K是常数)这样的线性表达式,然后就是分段使用矩阵快速密,最后得到结果。

接下来最重要的是怎么建立过渡矩阵:我们先来观察两个Fn,假设在i=3和4时P/n相等,设p=P/n,F3=C*A+D*B+p,F4=B*C+D*F3+p=B*C+C*D*A+D*D*B+D*p+p=(C+D*D)*B+(C*D)*A+(D+1)*p

则有F3=C*A+D*B+p   F4=(C+D*D)*B+(C*D)*A+(D+1)*p 那我们就要想办法将F3和F4的(A和B)的系数出现在同一个矩阵,首先我们先把A和B的初始系数确定好,设矩阵为a[3][3],a[0][0]=D(B的系数),a[0][1]=C(A的系数),因为要将D变成D*D+C,那么就要自乘一个D,再加上一个C,又因为B的系数是在第一行第一列,所以对第一行*第一列要得到D*D+C,那么a[1][0]=1,a[2][0]=0,而对于A的系数,则是第一行乘以第二列,要将C变成C*D,那么a[1][1]=a[2][1]=0,而对于p我们不妨设a[0][2]=p,而将P变成p*(D+1),那就是将第1行乘以第三列,那么a[1][2]=0,a[2][2]=1,到此矩阵就建立好了,而且不难发现在自乘之后,该矩阵的第二行就等于原矩阵的第一行,所以我们就能保证Fn-1和Fn-2的系数在同一个矩阵了。
#include<queue>
#include <cstdio>
#include<cstring>
#define debug(x) cout<<'x'<<' '<<x<<endl;
const int INF=0x3f3f3f3f;
typedef long long ll;
const int maxn = 1e5+10;
const int mod =1e9+7;
using namespace std;
int a,b,c,d,n,p,t;
struct mat
{
    int m[3][3];
    mat()
    {
        memset(m,0,sizeof(mat));
    }
    friend mat operator * (mat a,mat b)
    {
        mat c;
        for(int i=0;i<3;i++)
        {
            for(int j=0;j<3;j++)
            {
                ll t=0;
                for(int k=0;k<3;+k++)
                {
                    t+=(ll)a.m[i][k]*b.m[k][j];
                }
                c.m[i][j]=t%mod;
            }
        }
        return c;
    }
}I;
mat pow_mat(mat a,int b)
{
    mat c=I;
    while(b)
    {
        if(b&1){c=c*a;}
        a=a*a;
        b>>=1;
    }
    return c;
}
int main()
{
    I.m[0][0]=I.m[1][1]=I.m[2][2]=1;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d%d%d%d%d",&a,&b,&c,&d,&p,&n);
        if(n==1){printf("%d
",a);continue;}
        mat f;
        f.m[0][0]=d;
        f.m[0][1]=c;
        f.m[1][0]=1;
        f.m[2][2]=1;
        int f1=0;
        for(int i=3;i<=n;)
        {
            if(p/i==0)
            {
                mat w=f;
                w=pow_mat(w,n-i+1);
                ll ans=w.m[0][0]*(ll)b%mod+w.m[0][1]*(ll)a+w.m[0][2]%mod;
                ans%=mod;
                printf("%lld
",ans);
                f1=1;
                break;
            }
            int j=min(n,p/(p/i));
            mat w=f;
            w.m[0][2]=p/i;
            w=pow_mat(w,j-i+1);
            ll tmp1=(w.m[1][0]*(ll)b+w.m[1][1]*(ll)a+w.m[1][2])%mod;
            ll tmp2=(w.m[0][0]*(ll)b+w.m[0][1]*(ll)a+w.m[0][2])%mod;
            a=tmp1;b=tmp2;
            i=j+1;
        }
        if(!f1)printf("%d
",b);
    }
    return 0;
}
原文地址:https://www.cnblogs.com/kayiko/p/11284662.html