2021牛客多校 Day8 H.Scholomance Academy

题目描述

链接：https://ac.nowcoder.com/acm/contest/11259/H

Toilet-Ares finds a paper on the ground of Scholomance Academy, which says

$$ G(N)=\sum_{k_1+k_2+...+k_t=N}F(p_1^{k_1}p_2^{k_2}...p_t^{k_t}) $$

$$ F(n)=\sum_{a_1a_2...a_m=n}\varphi(a_1)\varphi(a_2)...\varphi(a_m) $$

​ Note that $\varphi(n)$ is the number of positive integer which is smaller than $n$ and co-prime with $n$. In particular $\varphi(1)=1$.

​ Now Toilet-Ares wonders, when given $N,t,m$ and $p_1,\ldots,p_t$, what the value of$G(N)$ modulo $998244353$ is supposed to be.

​ Of couse, Toilet-Ares who is phoning you, don't know how to solve this problem. Can you help him calculate the answer?

数据范围

$(1 \le N \le 10^9, 1 \le t \le 10^5, 1 \le m \le 5 , 1 \le p_i \le 10^5)$

It's guaranteed that $p_i$ is prime and different from each other

解题思路

​ 首先考虑$F(n)$ ，简单尝试后不难发现$F(n)$ 为积性函数，则$G(N)$ 可分解为 $\sum_{k_1+k_2+...+k_t=N}F(p_1^{k1})F(p_2^{k2})...F(p_t^{kt})$

简要证明：
设 $p,q$ 互质，则有

$$ F(p)F(q)=\sum_{a_1a_2...a_m=p}\varphi(a_1)\varphi(a_2)...\varphi(a_m)\sum_{b_1b_2...b_m=q}\varphi(b_1)\varphi(b_2)...\varphi(b_m) $$

$$ =\sum_{a_1a_2...a_mb_1b_2...b_m=pq}\varphi(a_1)\varphi(a_2)...\varphi(a_m)\varphi(b_1)\varphi(b_2)...\varphi(b_m) $$

令 $c_i=a_ib_i$ 则有 $\varphi(c_i)=\varphi(a_i)\varphi(b_i)$

$$ F(p)F(q)=\sum_{c_1c_2...c_m=pq}\varphi(c_1)\varphi(c_2)...\varphi(c_m)=F(pq) $$

​ 分解 $G(N)$ 后我们发现，对于 $F(n)$ ，我们只需求 $F(P_i^{k_i})$ 。 令 $f_{p_i}(x)=\sum_{k=0}^{\infty} F(p_i^{k})x^k$ 。则 $G(N)$ 就是 $G$ 的生成函数 $g(x)=\prod_{i=1}^t f_{p_i}(x)$ 的 $N$ 次项的系数。

​ 接下来考虑如何求出 $f_{p_i}(x)$ 。

​ 先考虑$f_{p_i}(x)$ 的 $k$ 次项系数 $F(p_i^k)$ 。由于 $p_i$ 为质数，于是我们可以改写 $F(p_i^k)$ 为

$$ F(p_i^k)=\sum_{d_1+d_2+...+d_m=k}\varphi(p_i^{d_1})\varphi(p_i^{d_2})...\varphi(p_i^{d_m}) $$

​ 再将其转换成生成函数的形式，$F(p_i^k)$ 即为 $h_{p_i}(x)=(\varphi(p_i^0)+\varphi(p_i^1)x+\varphi(p_i^2)x^2+...)^m$ 的 $k$ 次项的系数。 于是我们发现 $f_{p_i}(x)$ 就是 $h_{p_i}(x)$ ，并且根据欧拉函数的性质，当 $k\ge1$ 时，$\varphi(p_i^k)=(p_i-1)p_i^k$ 。

​ 代入得 $f_{p_i}(x)=(1+(p_i-1)p_ix+(p_i-1)p_i^2x^2+...)^m$

​ 发现这是一个等比数列，首项为 $(p_i-1)p_ix$ ，公比为 $p_ix$ 。

​ 使用求和公式，得到前 $n$ 项和为 $S_n=\frac{(p_i-1)p_i^{n+1}x^{n+1}-(p_i-1)x}{p_ix-1}$ 。这里令 $x\lt1$ ，当 $n\rightarrow\infty$ 时，$S=\frac{x-p_ix}{p_ix-1}$ 。

​ 于是

$$ f_{p_i}(x)=(1+\frac{x-p_ix}{p_ix-1})^m=\frac{(x-1)^m}{(p_ix-1)^m} $$

​ 直接代入 $g(x)$ ，得到

$$ g(x)=\prod_{i=1}^t\frac{(x-1)^m}{(p_ix-1)^m} $$

​ 更换顺序得

$$ g(x)=\frac{(\prod_{i=1}^tx-1)^m}{(\prod_{i=1}^tp_ix-1)^m} $$

​ 最后用二分递归求出 $\prod_{i=1}^tx-1$ 和 $\prod_{i=1}^tp_ix-1$ ，之后直接循环求出它们的 $m$ 次幂，再求这个分式的第 $N$ 次项的系数即可。

​ 求分式的第 $N$ 次项系数我使用波斯坦-茉莉算法，这个算法的优点是不需要任何矩阵知识，只使用到多项式乘法。复杂度 $O(k\log k \log n)$ ，$n$ 为多项式长度，$k$ 为要求的次项。

完整代码

#include<bits/stdc++.h>

using namespace std;
#define MAXN 100020
#define IOS std::ios::sync_with_stdio(false),std::cin.tie(0),std::cout.tie(0);

typedef long long ll;
typedef std::pair<int,int> pii;
typedef std::pair<long long,long long> pll;


const ll ALL_MOD=998244353;
namespace NFTS {
    const ll NFT_MOD = ALL_MOD, g = 3;
    std::vector<ll> rev, roots{0, 1};

    ll powMod(ll x, ll n) {
        ll res = 1;
        while (n) {
            if (n & 1) res = (res * x) % NFT_MOD;
            n >>= 1;
            x = (x * x) % NFT_MOD;
        }
        return res;
    }

    inline void DFT(std::vector<ll> &f) {
        int n = f.size();
        if (rev.size() != n) {
            int k = __builtin_ctz(n) - 1;
            rev.resize(n);
            for (int i = 0; i < n; i++) {
                rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
            }
        }
        if (roots.size() < n) {
            int k = __builtin_ctz(roots.size());
            roots.resize(n);
            while ((1 << k) < n) {
                ll e = powMod(g, ((NFT_MOD - 1) >> (k + 1)));
                for (int i = 1 << (k - 1); i < (1 << k); i++) {
                    roots[i << 1] = roots[i];
                    roots[i << 1 | 1] = (roots[i] * e) % NFT_MOD;
                }
                k++;
            }
        }
        for (int i = 0; i < n; i++) if (rev[i] < i) std::swap(f[i], f[rev[i]]);

        for (int k = 1; k < n; k <<= 1) {
            for (int i = 0; i < n; i += (k << 1)) {
                for (int j = 0; j < k; j++) {
                    ll u = f[i + j];
                    ll v = (f[i + j + k] * roots[k + j]) % NFT_MOD;
                    int x = u + v, y = u - v;
                    if (x >= NFT_MOD) x -= NFT_MOD;
                    if (y < 0) y += NFT_MOD;
                    f[i + j] = x;
                    f[i + j + k] = y;
                }
            }
        }
    }

    inline void IDFT(std::vector<ll> &f) {
        int n = f.size();
        std::reverse(f.begin() + 1, f.end());
        DFT(f);
        ll inv = powMod(n, NFT_MOD - 2);
        for (int i = 0; i < n; i++) {
            f[i] = (f[i] * inv) % NFT_MOD;
        }
    }
}// namespace NFTS


class Poly {
    void reverse() {
        std::reverse(data.begin(), data.end());
    }

public:
    static const ll POLY_MOD = ALL_MOD;
    static const ll INV2 = (POLY_MOD + 1) >> 1;
    std::vector<ll> data;

    Poly() {}

    Poly(int x) {
        if (x) data = {x};
    }

    Poly(const std::vector<ll> Data) : data(Data) {}

    Poly(const std::initializer_list<ll> &list) : data(list) {};

    inline const int size() const { return data.size(); };

    inline const ll operator[](const int &ind) const {
        return (ind < 0 || ind >= (int) data.size() ? 0ll : data[ind]);
    }

    inline ll &operator[](const int &ind) {
        return data[ind];
    };

    Poly operator-() const {
        auto b = data;
        for (auto &i:b) i = -i;
        return b;
    }

    Poly mulXn(const int &n) const {
        auto b = data;
        b.insert(b.begin(), n, 0);
        return b;
    }

    Poly modXn(const int &n) const {
        if (n > size()) return *this;
        return Poly({data.begin(), data.begin() + n});
    }

    Poly divXn(const int &n) const {
        if (size() <= n) return Poly();
        return Poly({data.begin() + n, data.end()});
    }

    Poly &operator+=(const Poly &b) {
        if (size() < b.size()) data.resize(b.size());
        for (int i = 0; i < b.size(); i++) {
            if ((data[i] += b.data[i]) > POLY_MOD) data[i] -= POLY_MOD;
        }
        return *this;
    }

    Poly &operator-=(const Poly &b) {
        if (size() < b.size()) data.resize(b.size());
        for (int i = 0; i < b.size(); i++) {
            if ((data[i] -= b.data[i]) < 0) data[i] += POLY_MOD;
        }
        return *this;
    }

    Poly &operator*=(Poly b) {
        int n = size(), m = b.size(), tot = std::max(1, m + n - 1);
        int sz = 1 << std::__lg((tot << 1) - 1);
        data.resize(sz);
        b.data.resize(sz);
        NFTS::DFT(data);
        NFTS::DFT(b.data);
        for (int i = 0; i < sz; i++) {
            data[i] = (data[i] * b.data[i]) % POLY_MOD;
        }
        NFTS::IDFT(data);
        data.resize(tot);
        return *this;
    }

    Poly &operator/=(Poly b) {
        int n = size(), m = b.size();
        if (n < m) return *this = Poly();
        reverse();
        b.reverse();
        *this *= b.inv(n - m + 1);
        data.resize(n - m + 1);
        reverse();
        return *this;
    }

    Poly &operator%=(const Poly &b) {
        return *this -= *this / b * b;
    }

    Poly operator+(const Poly &b) const {
        return Poly(*this) += b;
    }

    Poly operator-(const Poly &b) const {
        return Poly(*this) -= b;
    }

    Poly operator*(const Poly &b) const {
        return Poly(*this) *= b;
    }

    Poly operator/(const Poly &b) const {
        return Poly(*this) /= b;
    }

    Poly operator%(const Poly &b) const {
        return Poly(*this) %= b;
    }

    Poly powModPoly(int n, Poly p) const {
        Poly res(1), x(*this);
        while (n) {
            if (n & 1) (res *= x) %= p;
            n >>= 1;
            (x *= x) %= p;
        }
        return res;
    }

    ll inner(const Poly &b) {
        ll r = 0;
        int n = std::min(size(), b.size());
        for (int i = 0; i < n; i++) {
            r = (r + data[i] * b.data[i]) % POLY_MOD;
        }
        return r;
    }

    Poly derivation() const {
        if (data.empty()) return Poly();
        int n = size();
        std::vector<ll> r(n - 1);
        for (int i = 1; i < n; i++) r[i - 1] = (data[i] * i) % POLY_MOD;
        return r;
    }

    Poly integration() {
        if (data.empty()) return Poly();
        int n = size();
        std::vector<ll> r(n + 1), inv(n + 1, 1);
        for (int i = 2; i <= n; i++) inv[i] = (POLY_MOD - POLY_MOD / i) * inv[POLY_MOD % i] % POLY_MOD;
        for (int i = 0; i < n; i++) r[i + 1] = (data[i] * inv[i + 1]) % POLY_MOD;
        return r;
    }

    Poly inv(const int &n) const {
        assert(data[0] != 0);
        Poly x(NFTS::powMod(data[0], POLY_MOD - 2));
        int k = 1;
        while (k < n) {
            k <<= 1;
            x *= (Poly(2) - modXn(k) * x).modXn(k);
        }
        return x.modXn(n);
    }

    // 需保证首项为1
    Poly log(const int &n) const {
        return (derivation() * inv(n)).integration().modXn(n);
    }

    // 需保证首项为0
    Poly exp(const int &n) const {
        Poly x(1);
        int k = 1;
        while (k < n) {
            k <<= 1;
            x = (x * (Poly(1) - x.log(k) + modXn(k))).modXn(k);
        }
        return x.modXn(n);
    }

    // 需保证首项为1，开任意次方可log后exp
    Poly sqrt(const int &n) const {
        Poly x(1);
        int k = 1;
        while (k < n) {
            k <<= 1;
            x += modXn(k) * x.inv(k);
            x = x.modXn(k) * INV2;
        }
        return x.modXn(n);
    }

    // 减法卷积/转置卷积 {\rm MULT}(F(x),G(x))=\sum_{i\ge0}(\sum_{j\ge 0}f_{i+j}g_j)x^i
    Poly mulT(Poly b) const {
        if (b.size() == 0) return Poly();
        int n = b.size();
        b.reverse();
        return (*this * b).divXn(n - 1);
    }

    ll eval(const ll &x) const {
        ll r = 0, t = 1;
        int n = size();
        for (int i = 0; i < n; i++) {
            r = (r + data[i] * t) % POLY_MOD;
            t = (t * x) % POLY_MOD;
        }
        return r;
    }

    // 多点求值新科技：https://jkloverdcoi.github.io/2020/08/04/转置原理及其应用/
    std::vector<ll> evals(const std::vector<ll> &x) const {
        if (size() == 0) return std::vector<ll>(x.size());
        int n = x.size();
        std::vector<ll> ans(n);
        std::vector<Poly> g(4 * n);
        std::function<void(int, int, int)> build = [&](int l, int r, int p) {
            if (l == r) {
                g[p] = Poly({1, x[l] ? POLY_MOD - x[l] : 0});
            } else {
                int mid = (l + r) >> 1;
                build(l, mid, p << 1);
                build(mid + 1, r, p << 1 | 1);
                g[p] = g[p << 1] * g[p << 1 | 1];
            }
        };
        build(0, n, 1);

        std::function<void(int, int, int, const Poly &)> solve = [&](int l, int r, int p, const Poly &f) {
            if (r - l == 1) {
                ans[l] = f[0];
            } else {
                int mid = (l + r) >> 1;
                solve(l, mid, p << 1, f.mulT(g[p << 1 | 1]).modXn(mid - 1));
                solve(mid + 1, r, p << 1 | 1, f.mulT(g[p << 1]).modXn(mid - 1));
            }
        };
        solve(0, n, 1, mulT(g[1].inv(size())).modXn((n)));
        return ans;
    }

    // [x^k] f/g, bostan-mori O(k \log k \log n)  https://www.cnblogs.com/Potassium/p/15130342.html
    ll divat(Poly g, int k) const {
        Poly f = *this;
        int i;
        for (; k; k >>= 1) {
            Poly r = g;
            for (i = 1; i < r.size(); i += 2) r[i] = ALL_MOD - r[i];
            f *= r, g *= r;
            for (i = k & 1; i < f.size(); i += 2) f[i / 2] = f[i];
            f.data.resize(i / 2);

            for (i = 0; i < g.size(); i += 2) g[i / 2] = g[i];
            g.data.resize(i / 2);
        }
        return f.size() ? (f[0] * NFTS::powMod(g[0], ALL_MOD - 2)) % ALL_MOD : 0;

    }

    friend std::ostream &operator<<(std::ostream &out, const Poly &x) {
        for (auto &i:x.data) {
            out << i << ' ';
        }
        return out;
    }

};// class Poly


ll p[MAXN];

int n,t,m;

pair<Poly,Poly> solve(int l,int r){
    if(l==r){
        return {{ALL_MOD-1,1},{ALL_MOD-1,p[l]}};
    }
    int mid=(l+r)>>1;
    auto [a,b]=solve(l,mid);
    auto [c,d]=solve(mid+1,r);
    return {a*c,b*d};
}

int main(){
    IOS;
    
    cin>>n>>t>>m;

    for(int i=1;i<=t;i++) cin>>p[i];

    auto [f,g]=solve(1,t);
    Poly ff=f,gg=g;
    for(int i=1;i<m;i++){
        ff*=f,gg*=g;
    }

    cout<<ff.divat(gg,n);

    return 0;
}