P4389 付公主的背包
题意描述
付公主有一个可爱的背包qwq。
这个背包最多可以装 \(10^5\) 大小的东西
付公主有 \(n\) 种商品,她要准备出摊了
每种商品体积为 \(v_i\),都有无限件
给定 \(m\),对于 \(s\in [1,m]\),请你回答用这些商品恰好装 \(s\) 体积的方案数。
数据范围:\(n,m\leq 10^5,v_i\leq 10^5\)
solution
多项式+生成函数。
首先一件物品的生成函数显然为:\(\displaystyle 1+x^{v_i}+x^{2v_i}+x^{3v_i}+\cdots\) ,写成封闭形式为 \(\displaystyle 1\over 1-x^{v_i}\)
恰好装 \(s\) 体积的方案数的生成函数 \(F(x)\) ,显然是每个物品的生成函数卷起来之后的结果,即 \(\displaystyle F(x) = \prod_{i=1}^{n} {1\over 1-x^{v_i}}\)
考虑怎么求 \(F(x)\), 这玩意暴力展开的会比跑完全背包的复杂度还要差。
不妨考虑先求出 \(\ln F(x)\), 在 \(\exp\) 回去就可以得到 \(F(x)\) 。
显然有:\(\displaystyle \ln F(x) = \sum_{i=1}^{n} \ln({1\over 1-x^{v_i}})\) 。
那么我们的问题就转化为了怎么求 \(\displaystyle \ln({1\over1-x^{v_i}})\) 。
设 \(G(x) = \ln ({1\over 1-x^{v_i}})\), 对两边同时求导可得:
\(\displaystyle G^{\prime}(x) = \ln^{\prime}({1\over 1-x^{v_i}})\)
\(\displaystyle G^{\prime}(x) = (1-x^{v_i}){({1\over 1-x^{v_i}})^{\prime}}\)
我们把 \({1\over 1-x^{v_i}}\) 写成无穷项求和的形式来求导即:
\(\displaystyle ({1\over 1-x^{v_i}})^{\prime} = (\sum_{j=0}^{\infin} x^{jv_i})^{\prime} = \sum_{j=0}^{\infin} jv_i\times x^{jv_i-1}\)
把他带回去可得:
\(\displaystyle G^{\prime}(x) = \left(1-x^{v_i}\right) \sum_{j=0}^{\infin} jv_i\times x^{jv_i-1}\)
\(\displaystyle G^{\prime}(x) = \sum_{j=0}^{\infin} jv_i\times x^{jv_i-1} - \sum_{j=0}^{\infin} jv_i\times x^{(j+1)v_i-1}\)
因为 \(j = 0\) 的时候 \(jv_i\times x^{jv_i-1} = 0\) ,所以我们可以把第一个 \(sigma\) 中的下标改为从 \(1\) 开始,至于第二个 \(sigma\) 我们可以把下标平移,也变为了从 \(1\) 开始。
\(\displaystyle G^\prime(x) = \sum_{j=1}^{\infin} jv_i\times x^{jv_i-1} - \sum_{j=1}^{\infin} (j-1)v_i\times x^{jv_i-1}\)
\(\displaystyle G^{\prime}(x) = \sum_{j=1}^{\infin} v_i\times x^{jv_i-1}\)
把 \(G^{\prime}(x)\) 积分回去,就可以得到 \(G(x)\) 即:
\(\displaystyle G(x) = \int G^{\prime}(x)\,dx\)
\(\displaystyle G(x) = \int\left(\sum_{j=1}^{\infin} v_i\times x^{jv_i-1}\right)\,dx\)
\(\displaystyle G(x) =\sum_{j=1}^{\infin} {x^{jv_i}\over j}\)
则:\(\displaystyle \ln F(x) = \sum_{i=1}^{n}\sum_{j=0}^{\infin} {x^{jv_i}\over j}\)
对于每个 \(v_i\), 我们枚举它的倍数,给他倍数位置加上对应的系数,这样的复杂度是调和级数的。
求出来 \(\ln F(x)\) 之后,我们直接 \(\exp\) 回去即可。
Code
#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;
#define int long long
const int p = 998244353;
const int N = 1e6+10;
int n,m,x,rev[N],a[N],b[N],c[N],A[N],invA[N],invB[N],inv[N],v[N];
inline int read()
{
int s = 0,w = 1; char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-') w = -1; ch = getchar();}
while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
return s * w;
}
int ksm(int a,int b)
{
int res = 1;
for(; b; b >>= 1)
{
if(b & 1) res = res * a % p;
a = a * a % p;
}
return res;
}
void NTT(int *a,int lim,int opt)
{
for(int i = 0; i < lim; i++)
{
if(i < rev[i]) swap(a[i],a[rev[i]]);
}
for(int h = 1; h < lim; h <<= 1)
{
int wn = ksm(3,(p-1)/(h<<1));
if(opt == -1) wn = ksm(wn,p-2);
for(int j = 0; j < lim; j += (h<<1))
{
int w = 1;
for(int k = 0; k < h; k++)
{
int u = a[j + k];
int v = w * a[j + h + k] % p;
a[j + k] = (u + v) % p;
a[j + h + k] = (u - v + p) % p;
w = w * wn % p;
}
}
}
if(opt == -1)
{
int inv = ksm(lim,p-2);
for(int i = 0; i < lim; i++) a[i] = a[i] * inv % p;
}
}
void Inv(int n,int *a,int *b)
{
if(n == 1)
{
b[0] = ksm(a[0],p-2);
return;
}
Inv((n+1)>>1,a,b);
int lim = 1, tim = 0;
while(lim < (n<<1)) lim <<= 1, tim++;
for(int i = 0; i < lim; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(tim-1));
for(int i = 0; i < n; i++) c[i] = a[i];
for(int i = n; i < lim; i++) c[i] = 0;
NTT(b,lim,1); NTT(c,lim,1);
for(int i = 0; i < lim; i++) b[i] = (2 * b[i] % p - b[i] * b[i] % p * c[i] % p + p) % p;
NTT(b,lim,-1);
for(int i = n; i < lim; i++) b[i] = 0;
}
void qiudao(int n,int *a,int *b)
{
for(int i = 0; i < n-1; i++) b[i] = a[i+1] * (i+1) % p;
b[n-1] = 0;
}
void jifen(int n,int *a,int *b)
{
for(int i = 1; i < n; i++) b[i] = a[i-1] * inv[i] % p;
b[0] = 0;
}
void Ln(int n,int *a,int *b)
{
qiudao(n,a,A); Inv(n,a,invA);
int lim = 1, tim = 0;
while(lim < (n<<1)) lim <<= 1, tim++;
for(int i = 0; i < lim; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(tim-1));
NTT(A,lim,1); NTT(invA,lim,1);
for(int i = 0; i < lim; i++) A[i] = A[i] * invA[i] % p;
NTT(A,lim,-1); jifen(n,A,b);
for(int i = 0; i < lim; i++) A[i] = invA[i] = 0;
}
void Exp(int n,int *a,int *b)
{
if(n == 1)
{
b[0] = 1;
return;
}
Exp((n+1)>>1,a,b);
Ln(n,b,invB);
int lim = 1, tim = 0;
while(lim < (n<<1)) lim <<= 1, tim++;
for(int i = 0; i < lim; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(tim-1));
for(int i = 0; i < n; i++) c[i] = a[i];
for(int i = n; i < lim; i++) c[i] = 0;
NTT(invB,lim,1); NTT(b,lim,1); NTT(c,lim,1);
for(int i = 0; i < lim; i++) b[i] = (1 - invB[i] + c[i] + p) % p * b[i] % p;
NTT(b,lim,-1);
for(int i = n; i < lim; i++) b[i] = 0;
for(int i = 0; i < lim; i++) invB[i] = 0;
}
signed main()
{
n = read(); m = read();
for(int i = 1; i <= n; i++) v[read()]++;
for(int i = 1; i <= m; i++) inv[i] = ksm(i,p-2);
for(int i = 1; i <= m; i++)
{
if(!v[i]) continue;
for(int j = 1; i * j <= m; j++) a[i * j] = (a[i * j] + v[i] * inv[j] % p) % p;
}
Exp(m+1,a,b);
for(int i = 1; i <= m; i++) printf("%lld\n",b[i]);
return 0;
}