HDU - 6395 Sequence 矩阵快速幂 + 整除分块

传送门
除去常数,就是一个广义斐波那契数列,加上常数,就是把原来的二维矩阵变成三维矩阵
但是常数是在变化的,而利用整除分块的思想发现,在同一块的\(\frac{p}{i}\)是相同的。所以只需要维护两个值,\(fn_1\)表示该分块的第一个元素的前面第一个值,\(fn_2\)表示该分块的第一个元素的前面第二个值。
整除分块的左区间从3开始就行了,利用矩阵快速幂及时更新\(fn_1\)\(fn_2\)即可

#include <iostream>
#include <cstdio>
#include <cstring>
#define ll long long
using namespace std;
const int mod = 1e9 + 7;
const int N = 3;
struct Matrix{
    int n, m;
    ll a[N][N];
    Matrix(int n = 0, int m = 0):n(n),m(m){memset(a, 0, sizeof(a));};
    Matrix operator * (const Matrix &b) const {
        Matrix ans(n, b.m);
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < b.m; j++) {
                for(int k = 0; k < m; k++) {
                    ans.a[i][j] = (ans.a[i][j] + a[i][k] * b.a[k][j] % mod) % mod;
                }
            }
        }
        return ans;
    }
};
Matrix ksm(Matrix a, ll b){
    Matrix ans(a.n, a.m);
    for(int i = 0; i < max(a.n, a.m); i++) ans.a[i][i] = 1;
    while(b){
        if(b & 1) ans = ans * a;
        a = a * a;
        b >>= 1;
    }
    return ans;
}
ll a, b, c, d;
ll F(ll n, ll fn_1, ll fn_2, ll C_3) {
    Matrix base(3, 3);
    base.a[0][0] = d, base.a[0][1] = c, base.a[0][2] = 1;
    base.a[1][0] = 1;
    base.a[2][2] = 1;
    base = ksm(base, n - 2);
    Matrix ans(3, 1);
    ans.a[0][0] = fn_1;
    ans.a[1][0] = fn_2;
    ans.a[2][0] = C_3;
    ans = base * ans;
    return ans.a[0][0];
}
ll cal(ll n, ll p){ // fn_1 = fn - 1
    if(n == 1) return a;
    ll fn_1 = b, fn_2 = a, C_3 = p / 3;
    for(ll l = 3, r; l <= n; l = r + 1) {
        r = p / l ? min(p / (p / l), n) : n;
        C_3 = p / l;
        ll nowfn_1 = F(r - (l - 3), fn_1, fn_2, C_3);
        ll nowfn_2 = F(r - 1 - (l - 3), fn_1, fn_2, C_3);
        fn_1 = nowfn_1; fn_2 = nowfn_2;
    }
    return fn_1;
}
void solve(){
    ll n, p;
    scanf("%lld%lld%lld%lld", &a, &b, &c, &d);
    scanf("%lld%lld", &p, &n);
    printf("%lld\n", cal(n, p));
}  
int main(){
    int t;
    scanf("%d", &t);
    while(t--) solve();
    return 0;
}
posted @ 2020-10-26 17:03  Emcikem  阅读(91)  评论(0编辑  收藏  举报