【LOJ】#3090. 「BJOI2019」勘破神机

最新推荐文章于 2023-11-21 12:28:31 发布

转载最新推荐文章于 2023-11-21 12:28:31 发布 · 214 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/ivorysi/p/10984467.html

本文详细解析了洛谷平台上的题目“BJOI2019勘破神机”。通过学习BM算法，解决了斐波那契数列相关问题，利用特征根方程和多项式技巧，实现了高效求解。复杂度为k^2logr。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

LOJ#3090. 「BJOI2019」勘破神机

为了这题我去学习了一下BM算法。。

很容易发现这2的地方是$F_{1} = 1,F_{2} = 2$的斐波那契数列

3的地方是$G_{1} = 3,G_{2} = 11$其中下标表示长度的$\frac{1}{2}$，可以得到$G_{3} = 4G_{2} - G_{1}$

然后我们列一波特征根方程，可以得到

$m = 2$时
$$
\left{\begin{matrix}
x_{1} = \frac{1 + \sqrt{5}}{2}
\ x_{2} = \frac{1 - \sqrt{5}}{2}
\ c_{1} = \frac{5 + \sqrt{5}}{10}
\ c_{2} = \frac{5 - \sqrt{5}}{10}

\end{matrix}\right.
$$

$m = 3$时
$$
\left{\begin{matrix}
x_{1} = 2 + \sqrt{3}
\ x_{2} = 2 - \sqrt{3}
\ c_{1} = \frac{3 + \sqrt{3}}{6}
\ c_{2} = \frac{3 - \sqrt{3}}{10}

\end{matrix}\right.
$$
然后我们会可以把阶乘$n(n - 1)(n - 2)...(n - k + 1)$当成一个k次多项式

然后暴力算出每一项的系数，然后我们要算的就是

$\sum_{i = l}^{r} f_{i}^{k}$

然后把通项带进去

$\sum_{i = l}^{r} (c_{1}x_{1}^{i} + c_{2}x_{2}^{i})^{k}$

$\sum_{i = l}^{r}\sum_{j = 0}^{k}\binom{k}{j}c_{1}^{j}x_{1}^{ij}c_{2}^{k - j}x_{2}^{i(k - j)}$

$\sum_{j = 0}^{k}\binom{k}{j}c_{1}^{j}c_{2}^{k - j}\sum_{i = l}^{r}(x_{1}^{j}x_{2}^{k - j})^{i}$

后面是等比数列求和，算一下就好了

复杂度是$k^{2}\log r$的

#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int,int>
#define mp make_pair
#define pb push_back
#define space putchar(' ')
#define enter putchar('\n')
#define eps 1e-10
#define MAXN 2005
#define ba 47
//#define ivorysi
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
template<class T>
void read(T &res) {
    res = 0;T f = 1;char c = getchar();
    while(c < '0' || c > '9') {
    if(c == '-') f = -1;
    c = getchar();
    }
    while(c >= '0' && c <= '9') {
    res = res * 10 +c - '0';
    c = getchar();
    }
    res *= f;
}
template<class T>
void out(T x) {
    if(x < 0) {x = -x;putchar('-');}
    if(x >= 10) {
    out(x / 10);
    }
    putchar('0' + x % 10);
}
const int MOD = 998244353;
int inc(int a,int b) {
    return a + b >= MOD ? a + b - MOD : a + b;
}
int mul(int a,int b) {
    return 1LL * a * b % MOD;
}
void update(int &x,int y) {
    x = inc(x,y);
}
int fpow(int x,int c) {
    int res = 1,t = x;
    while(c) {
    if(c & 1) res = mul(res,t);
    t = mul(t,t);
    c >>= 1;
    }
    return res;
}
template<int T>
struct Complex {
    int r,i;
    Complex(int _r = 0,int _i = 0) {r = _r;i = _i;}
    friend Complex operator + (const Complex &a,const Complex &b) {
    return Complex(inc(a.r,b.r),inc(a.i,b.i));
    }
    friend Complex operator - (const Complex &a,const Complex &b) {
    return Complex(inc(a.r,MOD - b.r),inc(a.i,MOD - b.i));
    }
    friend Complex operator * (const Complex &a,const Complex &b) {
    return Complex(inc(mul(a.r,b.r),mul(T,mul(a.i,b.i))),inc(mul(a.r,b.i),mul(b.r,a.i)));
    }
    friend Complex fpow(Complex x,int64 c) {
    Complex res(1,0),t = x;
    while(c) {
        if(c & 1) res = res * t;
        t = t * t;
        c >>= 1;
    }
    return res;
    }
    friend Complex Query(Complex x,int64 n) {
    if(x.r == 1 && !x.i) return (n + 1)% MOD;
    Complex u = 1 - fpow(x,n + 1),d = 1 - x,t;
    t = d;t.i = MOD - t.i;
    int iv = fpow((t * d).r,MOD - 2);
    u = u * t * iv;
    return u;
    }
};
int m,k,f[1005],fac[1005],invfac[1005];
int64 l,r;
int C(int n,int m) {
    if(n < m) return 0;
    return mul(fac[n],mul(invfac[m],invfac[n - m]));
}
namespace mequal2 {
    Complex<5> x1,x2,c1,c2;
    void Main() {
    int inv2 = fpow(2,MOD - 2),inv10 = fpow(10,MOD - 2);
    x1 = Complex<5>(inv2,inv2);x2 = Complex<5>(inv2,MOD - inv2);
    c1 = Complex<5>(mul(5,inv10),inv10);c2 = Complex<5>(mul(5,inv10),MOD - inv10);
    int ans = mul(invfac[k],fpow((r - l + 1) % MOD,MOD - 2));
    Complex<5> res;
    for(int j = 1 ; j <= k ; ++j) {
        Complex<5> s;
        for(int h = 0 ; h <= j ; ++h) {
        Complex<5> t = C(j,h) * fpow(c1,h) * fpow(c2,j - h),q = fpow(x1,h) * fpow(x2,j - h),p;
        p = Query(q,r) - Query(q,l - 1);
        s = s + t * p;
        }
        res = res + s * f[j];
    }
    ans = mul(ans,res.r);
    out(ans);enter;
    }
}
namespace mequal3 {
    Complex<3> x1,x2,c1,c2;
    void Main() {
    int inv6 = fpow(6,MOD - 2);
    x1 = Complex<3>(2,1);x2 = Complex<3>(2,MOD - 1);
    c1 = Complex<3>(mul(3,inv6),inv6);c2 = Complex<3>(mul(3,inv6),MOD - inv6);
    int ans = mul(invfac[k],fpow((r - l + 1) % MOD,MOD - 2));
    r = r / 2;
    if(l & 1) l = (l + 1) / 2;
    else l = l / 2;
    
    Complex<3> res;
    if(l <= r) {
        for(int j = 1 ; j <= k ; ++j) {
        Complex<3> s;
        for(int h = 0 ; h <= j ; ++h) {
            Complex<3> t = C(j,h) * fpow(c1,h) * fpow(c2,j - h),q = fpow(x1,h) * fpow(x2,j - h),p;
            p = Query(q,r) - Query(q,l - 1);
            s = s + t * p;
        }
        res = res + s * f[j];
        }
    }
    ans = mul(ans,res.r);
    out(ans);enter;
    }
}
void Solve() {
    read(l);read(r);read(k);
    memset(f,0,sizeof(f));
    f[1] = 1;
    for(int i = 1 ; i < k ; ++i) {
    for(int j = i + 1 ; j >= 1 ; --j) {
        f[j] = inc(mul(MOD - i,f[j]),f[j - 1]);
    }
    }
    if(m == 2) mequal2::Main();
    else mequal3::Main();
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    int T;
    read(T);read(m);
    fac[0] = 1;
    for(int i = 1 ; i <= 1000 ; ++i) fac[i] = mul(fac[i - 1],i);
    invfac[1000] = fpow(fac[1000],MOD - 2);
    for(int i = 999 ; i >= 0 ; --i) invfac[i] = mul(invfac[i + 1],i + 1);
    for(int i = 1 ; i <= T ; ++i) {
    Solve();
    }
}

转载于:https://www.cnblogs.com/ivorysi/p/10984467.html