[HDU](6395)Sequence ---- 矩阵快速幂+分块

题目链接

做法：一开始推出了和正确题解一样的基础矩阵，可是自己想的太局限，没有想到分块的思想。还是缺乏练习呀！

具体做法，就是我们先把一部分的用O（1）的递推式计算好，然后后面的按照每相同的p/n 进行分块，然后这样就不会错了。

 参考题解：https://www.cnblogs.com/DynastySun/p/9470717.html

#include<bits/stdc++.h>
#define rep(i,s,t) for(int i = (int)(s); i <= (int)(t); i++)
#define rev(i,t,s) for(int i = (int)(t); i >= (int)(s); i--)
#define pb(x) push_back(x)
#define all(x) x.begin(),x.end()
#define sz(x) (int)(x).size()
using namespace std;
typedef long long ll;
const int mod = 1e9+7;
const double PI = 4*atan(1.0);
const int maxn = 5e5+5;
const int INF = 0x3f3f3f3f;
ll A,B,C,D,P,N;
int f[maxn];
struct mat{
    ll a[15][15];
    mat(){
        memset(a,0,sizeof(a));
    }
};
mat mul(mat x,mat y)
{
    mat res;
    for(int i=0;i<3;i++)
        for(int j=0;j<3;j++)
            for(int k=0;k<3;k++)
                res.a[i][j] = (res.a[i][j]+x.a[i][k]*y.a[k][j])%mod;
    return res;
}
ll qpow(ll fn_1, ll fn_2, ll c, ll b)
{
    mat res;
    for(int i=0;i<3;i++)
        res.a[i][i] = 1;
    mat bas;
    bas.a[0][0] = D;bas.a[0][1] = C;bas.a[0][2] = 1;
    bas.a[1][0] = 1;bas.a[1][1] = 0;bas.a[1][2] = 0;
    bas.a[2][0] = 0;bas.a[2][1] = 0;bas.a[2][2] = 1;
    while(b)
    {
        if(1&b) res = mul(res,bas);
        bas = mul(bas,bas);
        b>>=1;
    }
    ll ans = (res.a[0][0]*fn_1%mod+res.a[0][1]*fn_2%mod+res.a[0][2]*c%mod)%mod;
    return ans;
}
void Init()
{
    f[1] = A;
    f[2] = B;
    for(int i=3;i<maxn;i++)
        f[i] = (C*f[i-2]%mod+D*f[i-1]%mod+(P/i))%mod;
}
int main()
{
    #ifdef LOCAL_FILE
    freopen("in.txt","r",stdin);
    #endif // LOCAL_FILE
    ios_base::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    int t;
    cin>>t;
    while(t--)
    {
        ll ans = 0;
        cin>>A>>B>>C>>D>>P>>N;
        Init();
        if(N<maxn)
            cout<<f[N]<<endl;
        else
        {
            ll m = maxn;
            ll a = f[m-1];
            ll b = f[m-2];
            ll c = P/m;
            ll num;
            if(c!=0)//从第m项开始分块计算
            {
                num = P/c-m+1;//计算出P/m = c的有多少项
                while(m+num-1<N)//计算出f[m]~f[n]项
                {
                    ll tfn_,tfn_1;
                    tfn_ = qpow(a,b,c,num);//算出中间的fm,num个相同的[p/m]
                    tfn_1 = qpow(a,b,c,num-1);
                    m+=num;
                    a = tfn_;//更新f[m-1]
                    b = tfn_1;//更新f[m-2]
                    c = P/m;
                    if(c == 0) break;//c = 0 说明后面的c都为0了，剩下的直接用快速幂
                    num = P/c-m+1;
                }
            }
            num = N-m+1;
            ans = qpow(a,b,c,num);
            cout<<ans<<endl;
        }
    }
    return 0;
}