【2021牛客暑期多校训练营8H】Scholomance Academy

题面蛮长的，大概是定义了两个函数

嗯，t(1e5)和m(1-5)还有那些个素数是给定的，现在给你个1e9的N，让你算。
那么第一很明显的是，欧拉函数是积性函数，而的公式是个迪利克雷卷积形式，所以显然也是积性函数，所以我们就可以得到

然后这个形式他就很想让人用生成函数。根据欧拉函数的性质.

因此

其实到这里都没什么难度，我们只要求出的系数就好了。题解说是线性递推，但是很多人（好吧只有我）不知道为什么能转化成线性递推。

然后那天在群里问，有人说展开，我会想起了斐波那契的生成函数解法，于是才想到。
事实上是这样的将生成函数上线展开，设成这样。

然后移过来

(其实到这里参考生成函数求斐波那契，大家应该会了)也就是

我们假设

考虑等式（1）两边的每项系数

...

...

到这里你应该看出来了，的前项系数可以用多项式求逆（对）和的卷积得到，而后面可以用线性递推得到。
后面就不写了，大家可以参考很多地方的齐次线性递推的blog or 板子，就很简单。

最后放一发队友的代码（多项式求逆，取模，递推，快速幂板子我全是抄的逆十字的所以就不放了）

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef vector<int> vi;
#define endl '\n'
// const int mod = 1e9 + 7;
const int mod = 998244353, G = 3, Gi = 332748118;
const int inf = 1e9;
const double pi = acos(-1);
#define int ll

typedef vector<int> poly;

int lim;
int rev[2100000], p[100005], ip[500005];
int n, t, m;

void invSieve(int n, int p) {
    ip[0] = ip[1] = 1;
    for (int i = 2; i <= n; i++) {
        ip[i] = (ll) (p - p / i) * ip[p % i] % p;
    }
}

ll ksm(ll a, ll x) {
    ll ans = 1;
    while (x) {
        if (x & 1) ans = ans * a % mod;
        a = a * a % mod;
        x /= 2;
    }
    return ans;   
}

ll inv(int x) {
    return ksm(x, mod - 2);
}

void init(int n) {
    for(lim = 1; lim < n; lim <<= 1);
    for(int i = 0; i < lim; ++i)
        rev[i] = rev[i >> 1] >> 1 | ((i & 1) ? lim >> 1 : 0);
}

int sub(int a, int b) {return a < b ? a - b + mod : a - b; }

// TODO: NTT
void NTT(int n, poly &a, int t) {
    for (int i = 0; i < n; ++i) {
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    }
    for (int len = 1; len <= (n >> 1); len <<= 1) {
        int wn = ksm(t == 1 ? G : Gi, (mod - 1) / (len * 2));
        for (int i = 0; i <= n - (len << 1); i += (len << 1)) {
            int w = 1;
            for (int j = 0; j < len; ++j) {
                int x = a[i + j], y = 1ll * w * a[i + j + len] % mod;
                a[i + j] = (x + y) % mod;
                a[i + j + len] = sub(x, y);
                w = 1ll * w * wn % mod;
            }
        }
    }
    if (t == -1) {
        int ki = inv(n);
        for(int i = 0; i < n; ++i) a[i] = 1ll * a[i] * ki % mod;
    }
}

poly operator * (poly a, int b) {
    for (int i = 0; i < a.size(); ++i) a[i] = 1ll * a[i] * b % mod;
    return a;
}

poly operator + (poly a, poly b) {
    poly c = a;
    c.resize(max(a.size(), b.size()));
    for (int i = 0; i < b.size(); ++i) c[i] = (c[i] + b[i]) % mod;
    return c;
}

poly operator - (poly a, poly b) {
    poly c = a;
    c.resize(max(a.size(), b.size()));
    for (int i = 0; i < b.size(); ++i) c[i] = sub(c[i], b[i]);
    return c;
}

// TODO: 多项式乘法
poly operator * (poly a, poly b) {
    int deg = a.size() + b.size() - 1;
    if (a.size() <= 32 || b.size() <= 32) {
        poly c(deg, 0);
        for (int i = 0; i < a.size(); ++i)
            for (int j = 0; j < b.size(); ++j)
                (c[i + j] += 1ll * a[i] * b[j] % mod) %= mod;
        return c;
    }
    init(deg);
    a.resize(lim), b.resize(lim);
    NTT(lim, a, 1), NTT(lim, b, 1);
    for (int i = 0; i < lim; ++i) a[i] = 1ll * a[i] * b[i] % mod;
    NTT(lim, a, -1);
    a.resize(deg);
    return a;
}

// TODO: 多项式微分
poly Dif(poly a) {
    for (int i = 0; i + 1 < a.size(); i++) a[i] = 1ll * a[i + 1] * (i + 1) % mod;
    a.pop_back();
    return a;
}

// TODO: 多项式积分
poly Int(poly a) {
    for (int i = a.size() - 1; i > 0; i--) a[i] = 1ll * a[i - 1] * inv(i) % mod;
    a[0] = 0;
    return a;
}

// TODO: 多项式求逆
poly Inv(poly a, int deg = -1) {
    if (deg == -1) deg = a.size();
    poly b(1, inv(a[0])), c;
    for (int len = 2; (len >> 1) < deg; len <<= 1) {
        c.resize(len);
        init(len << 1);
        for (int i = 0; i < len; i++) c[i] = i < a.size() ? a[i] : 0;
        c.resize(lim), b.resize(lim);
        NTT(lim, c, 1), NTT(lim, b, 1);
        for (int i = 0; i < lim; i++) b[i] = b[i] * (998244355 - 1ll * c[i] * b[i] % mod) % mod;
        NTT(lim, b, -1);
        b.resize(len);
    }
    b.resize(deg);
    return b;
}

// TODO: 多项式ln
poly Ln(poly a, int deg = -1) {
    if (deg == -1) deg = a.size();
    a = Dif(a) * Inv(a, deg);
    a.resize(deg);
    a = Int(a);
    return a;
}

// TODO: 多项式exp
poly Exp(poly a, int deg = -1) {
    if (deg == -1) deg = a.size();
    poly b(1, 1), c;
    int n = a.size();
    for (int len = 2; (len >> 1) < deg; len <<= 1) {
        c = Ln(b, len);
        for (int i = 0; i < len; i++) {
            c[i] = sub((i < n ? a[i] : 0), c[i]);
        }
        c[0]++;
        b = b * c;
        b.resize(len);
    }
    b.resize(deg);
    return b;
}

// TODO: 多项式除法
poly operator / (poly a, poly b) {
    int deg = a.size() - b.size() + 1;
    reverse(a.begin(), a.end());
    a.resize(deg);
    reverse(b.begin(), b.end());
    a = a * Inv(b, deg);
    a.resize(deg);
    reverse(a.begin(), a.end());
    return a;
}

// TODO: 多项式取模
poly operator % (poly a, poly b) {
    if (a.size() < b.size()) return a;
    a = a - (b * (a / b));
    a.resize(b.size() - 1);
    return a;
}

poly solve(int l, int r) {
    if (l == r) {
        poly a = {1, mod - p[(l - 1) / m + 1]};
        return a;
    }
    int mid = (l + r) / 2;
    return solve(l, mid) * solve(mid + 1, r);
}

signed main() {
    ios::sync_with_stdio(false); cin.tie(0);
#ifndef ONLINE_JUDGE
    freopen("a.in", "r", stdin);
    freopen("a.out", "w", stdout);
#endif
    cin >> n >> t >> m;
    int N = m * t;
    poly a;
    int C = 1, f = 1;
    for (int i = 1; i <= t; i++) cin >> p[i];
    invSieve(N, mod);
    for (int i = 0; i <= N; i++, f *= -1) {
        a.push_back(f == 1 ? C : mod - C);
        C = 1ll * C * (N - i) % mod * ip[i + 1] % mod;
    }
//     for (int i : a) cout << i << ' ';
//     cout << endl;
    poly b = solve(1, N);
//     for (int i : b) cout << i << ' ';
//     cout << endl;
    poly c = a * Inv(b);
//     for (int i : c) cout << i << ' ';
//     cout << endl;
    c.resize(N + 1);
    if (n <= N) cout << (c[n] + mod) % mod << endl, exit(0);
    reverse(b.begin(), b.end());
    n--;
    poly t0 = {1}, t1 = {0, 1};
    while (n) {
        if (n & 1) t0 = t0 * t1 % b;
        t1 = t1 * t1 % b;
        n /= 2;
    }
    int ans = 0;
    t0.resize(N);
    for (int i = 0; i < N; i++) {
        ans = (ans + 1ll * c[i + 1] * t0[i]) % mod;
    }
    cout << (ans + mod) % mod << endl;
}