luogu4389：付公主的背包

luogu4389：付公主的背包

题意：

给\(n\leq 10^5\)个东西，每个物品有体积\(v_i\)，有无限件。

给定\(m\)，对于\(s\in [1,m]\)，问用这些物品恰好装\(s\)体积的方案数。

对\(998244353\)取模。

思路：

如果数据是\(10^4\)，那么这就是一道完全背包的题目。

但显然此时我们需要效率更高的算法。

我们知道在求凑成某种数值的方案数的时候可以先设计母函数。

\[G(x)=\prod(1+x^{v_i}+x^{2v_i}+...) \]

这个公式展开后，对于\(s\in[1,m]\)的项数的系数就是对应方案数。

暴力展开的话，时间复杂度会达到\(O(n^2logn)\)，还不如写完全背包。

考虑优化。

第\(i\)个物品的多项式为：

\[f(i)=\sum_{j=0}x^{jv_i}=\frac{1}{1-x^{v_i}} \]

可能有学过收敛域的同学可能会觉得有些疑惑，\(p\)级数要求公比的绝对值严格小于\(1\)，这里其实\(x\)只是一个符号而已，取值自定，所以我可以假设他是小于\(1\)的公比。

所以我们需要的多项式就是：

\[g=\prod_{i=1}^nf(i) \]

改乘法为加法：

\[ln(g)=\prod ln(f(i))=ln(\sum f(i)) \]

那么多项式如何求\(ln\)呢？

\[lnf=\int \frac{f'}{f} \]

我们把\(f(i)\)带入，下面用\(v\)代替\(v_i\)。

\[lnf(i)=ln\frac{1}{1-x^v}=\int \frac{vx^{v-1}}{1-x^v}dx \]

将\(f(i)=\frac{1}{1-x^v}\)代入：

\[\int vx^{v-1}f(i)=\int\sum_{i=0}^\infin x^{iv}\times vx^{v-1}=\int\sum_{i=1}^\infin vx^{iv-1} \]

积分：

\[\sum_{i=1}^\infin \frac{1}{i}x^{iv} \]

这也就是\(ln(\sum f(i))\)了。

但是我们要求的是\(\prod f(i)\)，所以还需要\(exp\)回来。

代码如下：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int G = 3;
const int mod = 998244353;
const int maxn = 5e5+10;

ll qmi(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if(b&1) res = (res*a)%mod;
        b >>= 1;
        a = (a*a)%mod;
    } return res%mod;
}

int limit, len, rev[maxn];
void NTT(ll *c, int op)
{
    for(int i = 1; i <= limit; i++)
        if(i < rev[i]) swap(c[i], c[rev[i]]);

    for(int mid = 1; mid < limit; mid <<= 1)
    {
        ll gn = qmi(G, (mod-1)/(mid<<1));
        if(op == -1) gn = qmi(gn, mod-2);
        for(int j = 0, R = mid<<1; j < limit; j += R)
        {
            ll g = 1;
            for(int k = 0; k < mid; k++, g = (g*gn)%mod)
            {
                ll x = c[j+k], y = g*c[j+k+mid]%mod;
                c[j+k] = (x+y)%mod;
                c[j+k+mid] = (x-y+mod)%mod;
            }
        }
    }
    if(op == 1) return;
    ll inv = qmi(limit, mod-2);
    for(int i = 0; i <= limit; i++) c[i] = c[i]*inv%mod;
}

ll tmp[maxn];
void get_inv(ll *a, ll *b, int deg)
{
	if(deg == 1)
	{
		b[0] = qmi(a[0], mod-2);
		return;
	}
	get_inv(a, b, (deg+1)>>1);

	len = 0, limit = 1;
	while(limit <= (deg+deg)) limit <<= 1, len++;

	for(int i = 1; i <= limit; i++)
		rev[i] = (rev[i>>1]>>1)|((i&1)<<(len-1));

	for(int i = 0; i < deg; i++) tmp[i] = a[i];
	for(int i = deg; i <= limit; i++) tmp[i] = 0;

	NTT(tmp, 1), NTT(b, 1);
	for(int i = 0; i <= limit; i++)
		b[i] = 1ll*(2-1ll*tmp[i]*b[i]%mod+mod)%mod*b[i]%mod;
	NTT(b, -1);
	for(int i = deg; i <= limit; i++) b[i] = 0;
}

void deriv(ll *a, int n)
{
    for(int i = 1; i < n; i++)
        a[i-1] = i*a[i]%mod;
    a[n-1] = 0;
}

void integ(ll *a, int n)
{
    for(int i = n-1; i > 0; i--)
        a[i] = a[i-1]*qmi(i, mod-2)%mod;
    a[0] = 0;
}

void get_ln(ll *a, ll *b, int deg)
{
    get_inv(a, b, deg);
    for(int i = 0; i < deg; i++)
        tmp[i] = a[i];
    for(int i = deg; i <= limit; i++)
        tmp[i] = 0;
    deriv(tmp, deg);
    NTT(tmp, 1); NTT(b, 1);
    for(int i = 0; i <= limit; i++)
        b[i] = b[i]*tmp[i]%mod;
    NTT(b, -1);
    integ(b, deg);
}

ll c[maxn];
void get_exp(ll *a, ll *b, int n)
{
    if(n == 1) {
        b[0] = 1;
        return;
    }
    get_exp(a, b, (n+1)>>1);
    memset(c, 0, sizeof c);
    get_ln(b, c, n); // c = lnb
    c[0] = (1-c[0]+a[0]+mod)%mod;
    for(int i = 1; i < n; i++)
        c[i] = (a[i]-c[i]+mod)%mod;
    NTT(c, 1); NTT(b, 1);
    for(int i = 0; i <= limit; i++)
        b[i] = b[i]*c[i]%mod;
    NTT(b, -1);
    for(int i = n; i <= limit; i++)
        b[i] = 0;
}

int n, m;
int cnt[maxn];
ll f[maxn], g[maxn];

int main()
{
    scanf("%d%d", &n, &m);
    for(int i = 1, x; i <= n; i++){
        scanf("%d", &x);
        cnt[x]++;
    }

    for(int v = 1; v <= m; v++)
    if(cnt[v])
    {
        for(int iv= v; iv <= m; iv += v)
        f[iv]=(f[iv]+cnt[v]*qmi(iv/v,mod-2)%mod)%mod;
    }

    get_exp(f, g, m+1);
    for(int i = 1; i <= m; i++)
        printf("%lld\n", g[i]);
    return 0;
}