洛谷P4389 付公主的背包

设商品体积为的为，

则所求答案为

麦克劳林公式展开后可得

枚举倍数后可得系数，再求一次多项式即可得到答案，时间复杂度

#include <bits/stdc++.h>
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/trie_policy.hpp>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/priority_queue.hpp>

using namespace std;
using namespace __gnu_pbds;
using namespace __gnu_cxx;
typedef long long ll;
//tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> tr;
//__gnu_pbds::priority_queue<int, greater<int>, pairing_heap_tag> qu;
//typedef trie<string, null_type, trie_string_access_traits<>, pat_trie_tag, trie_prefix_search_node_update> pref_trie;
const int N = 200005;
int n, m, k;
namespace polybase {//范围为1e9需要先取模,别忘记改模数
    const ll mod = 998244353;
    int limit;
    ll _wn[25];

    ll fpow(ll x, ll r)
    {
        ll result = 1;
        while (r)
        {
            if (r & 1)result = result * x % mod;
            r >>= 1;
            x = x * x % mod;
        }
        return result;
    }

    int _ = []
    {
        for (int i = 0; i <= 23; i++)_wn[i] = fpow(3, (mod - 1) >> i);
        return 0;
    }();

    inline int norm(int n) { return 1 << __lg(n * 2 - 1); }

    void NTT(ll *A, int type)
    {
        int i, j = limit >> 1, k, l, c = 0;
        ll u, v, w, wn;
        for (i = 1; i < limit - 1; i++)
        {
            if (i < j)swap(A[i], A[j]);
            for (k = limit >> 1; (j ^= k) < k; k >>= 1);
        }

        for (l = 2; l <= limit; l <<= 1)
        {
            i = l >> 1, wn = _wn[++c];
            for (j = 0; j < limit; j += l)
            {
                w = 1;
                for (k = j; k < j + i; k++)
                {
                    u = A[k], v = A[k + i] * w % mod;
                    A[k] = u + v;
                    if (A[k] >= mod)A[k] -= mod;

                    A[k + i] = u - v;
                    if (A[k + i] < 0)A[k + i] += mod;
                    w = w * wn % mod;
                }
            }
        }
        if (type == -1)
        {
            ll inv = fpow(limit, mod - 2);
            for (i = 0; i < limit; i++)A[i] = A[i] * inv % mod;
            for (i = 1; i < limit / 2; i++)swap(A[i], A[limit - i]);
        }
    }

    struct poly : public vector<ll>
    {
        using vector<ll>::vector;
#define T (*this)

        poly modxk(int k) const
        {
            k = min(k, (int) size());
            return poly(begin(), begin() + k);
        }

        poly rev() const { return poly(rbegin(), rend()); }

        friend void NTT(poly &a, const int type) { NTT(a.data(), type); }

        friend poly operator*(const poly &x, const poly &y)
        {
            if (x.empty() || y.empty())return poly();
            poly a(x), b(y);
            int len = a.size() + b.size() - 1;
            limit = norm(len);
            a.resize(limit), b.resize(limit);
            NTT(a, 1);
            NTT(b, 1);
            for (int i = 0; i < limit; i++)
                a[i] = a[i] * b[i] % mod;
            NTT(a, -1);
            a.resize(len);
            return a;
        }

        poly operator+(const poly &b)
        {
            poly a(T);
            if (a.size() < b.size())
                a.resize(b.size());
            for (int i = 0; i < b.size(); i++)//用b.size()防止越界
            {
                a[i] += b[i];
                if (a[i] >= mod)a[i] -= mod;
            }
            return a;
        }

        poly operator-(const poly &b)
        {
            poly a(T);
            if (a.size() < b.size())
                a.resize(b.size());
            for (int i = 0; i < b.size(); i++)
            {
                a[i] -= b[i];
                if (a[i] < 0)a[i] += mod;
            }
            return a;
        }

        poly operator*(const ll p)
        {
            poly a(T);
            for (auto &x:a)
                x = x * p % mod;
            return a;
        }

        poly &operator<<=(int r) { return insert(begin(), r, 0), T; }//注意逗号,F(x)*(x^r)

        poly operator<<(int r) const { return poly(T) <<= r; }

        poly operator>>(int r) const { return r >= size() ? poly() : poly(begin() + r, end()); }

        poly &operator>>=(int r) { return T = T >> r; }//F[x]/(x^r)

        poly deriv()//求导
        {
            if (empty())return T;
            poly a(size() - 1);
            for (int i = 1; i < size(); i++)//注意是size()
                a[i - 1] = T[i] * i % mod;
            return a;
        }

        poly integ()//积分
        {
            poly a(size() + 1);
            for (int i = 1; i < a.size(); i++)//注意是a.size()
                a[i] = T[i - 1] * fpow(i, mod - 2) % mod;
            return a;
        }

        poly inv(int n)
        {
            poly a{fpow(T[0], mod - 2)};
            int k = 1;
            while (k < n)
            {
                k <<= 1;
                a = (a * 2 - modxk(k) * a * a).modxk(k);
            }
            return a.modxk(n);
        }

        poly sqrt(int n)//f[0]必须等于1
        {
            poly a{1};
            int k = 1;
            const ll inv2 = fpow(2, mod - 2);
            while (k < n)
            {
                k <<= 1;
                a = ((modxk(k) * a.inv(k)).modxk(k) + a) * inv2;
            }
            return a.modxk(n);
        }

        poly ln(int n)//需要保证f[0]=1
        {
            return (deriv() * inv(n)).integ().modxk(n);
        }

        poly exp(int n)//需要保证f[0]=1
        {
            poly a{1};
            int k = 1;
            while (k < n)
            {
                k <<= 1;
                a = (a * (poly{1} - a.ln(k) + modxk(k))).modxk(k);
            }
            return a.modxk(n);
        }

#undef T
    };
}
using namespace polybase;
int cnt[N];
ll inv[N];

int main()
{
    int p, q, u, v, w, x, y, z, T;
    cin >> n >> m;
    for (int i = 1; i <= n; i++)
        scanf("%d", &x), cnt[x]++;
    for (int i = 1; i <= m; i++)
        inv[i] = fpow(i, mod - 2);
    poly F(m + 1);
    for (int i = 1; i <= m; i++)
        for (int j = 1; i * j <= m; j++)
            F[i * j] = (F[i * j] + inv[j] * cnt[i]) % mod;
    F = F.exp(m + 1);
    for (int i = 1; i <= m; i++)
        printf("%d\n", F[i]);
    return 0;
}