2018 多校 HDU

这题就是一个矩阵快速幂，但是后面还有一项，重点在对这个项的处理上面，关于这个项的处理，可以看另一篇blog（https://www.cnblogs.com/DynastySun/p/9462847.html），这上面有关于那个项的一些特点。我们可以观察多，这个项的变化规律，P/n（向下取整）,当n越小的时候，这个值变化的越频繁，当n越大，这个值变化的越小。这就好办了，对于变化频繁的一段我们直接使用递推式子，将前50000项打表打出来，后面的变化不那么频繁的，我们就使用矩阵快速幂求解就可以了。这个题就是代码会有点多（重要部分，代码中已经给出解释）：

#include <bits/stdc++.h>
using namespace std;

typedef long long LL;

const LL MOD = 1e9+7;
const int maxn = 500000;

LL A, B, C, D, n, P, F[maxn];

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;
}

void mut(LL a[][3], LL b[][3]){
    LL t[3][3];
    memset(t, 0, sizeof(t));
    for(int i = 0; i < 3; ++i)
        for(int j = 0; j < 3; ++j)
            for(int k = 0; k < 3; ++k)
                t[i][j] = (t[i][j] + a[i][k]*b[k][j]%MOD)%MOD;
    memcpy(a, t, sizeof(t));
}

LL quick_pow(LL a, LL b, LL c, LL m)
{
    LL ans[][3] = {{1,0,0},{0,1,0},{0,0,1}};
    LL t[][3] = {{D,C,1},{1,0,0},{0,0,1}};
    while(m){
        if(m&1) mut(ans, t);
        mut(t, t);
        m >>= 1;
    }
    return (ans[0][0]*a%MOD + ans[0][1]*b%MOD + ans[0][2]*c)%MOD;
}

int main(){
    int T;
    cin >> T;
    while(T--)
    {
        cin >> A >> B >> C >> D >> P >> n;
        init();
        if(n < maxn)
        {
            printf("%lld
", F[n]);
        }
        else
        {
            LL m = maxn;    // 计算 f[m]
            LL c = P/m;     // 计算 P/m
            LL a = F[m-1];  // 存下 f[m-1] 用于矩阵快速幂
            LL b = F[m-2];  // 存下 f[m-2] 用于矩阵快速幂
            LL num;         // 代表 P/n == c 的n的数量
            if(c != 0)      // 如果 c == 0 我们就可以直接计算了
            {
                num = P/c - m + 1;  // 计算这之间的数量
                // 我们本次计算可以算到f[m+num-1]，这就要求n+num-1 < n
                // 我们要最后一次矩阵快速幂刚好能直接算到第n项
                while(m + num <= n)
                {
                    LL ta, tb;
                    ta = quick_pow(a, b, c, num);
                    if(num == 1)
                        tb = a;
                    else
                        tb = quick_pow(a, b, c, num-1);
                    m += num;
                    a = ta;         //  更新目前的f[m-1]
                    b = tb;         //  更新目前的f[m-2]
                    c = P/m;
                    if(c == 0) break;
                    num = P/c - m + 1;
                }
            }
            num = n - m + 1;
            LL ans = quick_pow(a, b, c, num);
            printf("%lld
", ans);
        }
    }
    return 0;
}