Luogu4389 付公主的背包(生成函数+多项式exp)
显然构造出生成函数,对体积v的物品,生成函数为1+xv+x2v+……=1/(1-xv)。将所有生成函数乘起来得到的多项式即为答案,设为F(x),即F(x)=1/∏(1-xvi)。但这个多项式的项数是Σvi级别的,无法直接分治FFT卷起来。
我们要降低多项式的次数,于是考虑取对数,化乘为加,得到lnF(x)=-Σln(1-xvi)。只要对每个多项式求出ln加起来再exp回去即可。
考虑怎么对这个特殊形式的多项式求ln。对ln(1-xv)求导,得ln(1-xv)'=(1-xv)'/(1-xv)=-vxv-1/(1-xv)=-vΣx(k+1)v-1,再积分得-vΣx(k+1)v/(k+1)v=-Σxkv/k(k>=1)。注意到由调和级数,总共只有mlogm项要加起来。然后多项式exp即可得到F(x)。
#include<iostream> #include<cstdio> #include<cmath> #include<cstdlib> #include<cstring> #include<algorithm> using namespace std; #define ll long long #define P 998244353 #define N 140010 char getc(){char c=getchar();while ((c<'A'||c>'Z')&&(c<'a'||c>'z')&&(c<'0'||c>'9')) c=getchar();return c;} int gcd(int n,int m){return m==0?n:gcd(m,n%m);} int read() { int x=0,f=1;char c=getchar(); while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();} while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar(); return x*f; } int n,m,t,a[N],f[N<<2],g[N<<2],x[N<<2],y[N<<2],r[N<<2],A[N<<2],B[N<<2]; int ksm(int a,int k) { int s=1; for (;k;k>>=1,a=1ll*a*a%P) if (k&1) s=1ll*s*a%P; return s; } int inv(int a){return ksm(a,P-2);} void DFT(int *a,int n,int g) { for (int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|(i&1)*(n>>1); for (int i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]); for (int i=2;i<=n;i<<=1) { int wn=ksm(g,(P-1)/i); for (int j=0;j<n;j+=i) { int w=1; for (int k=j;k<j+(i>>1);k++,w=1ll*w*wn%P) { int x=a[k],y=1ll*w*a[k+(i>>1)]%P; a[k]=(x+y)%P,a[k+(i>>1)]=(x-y+P)%P; } } } } void IDFT(int *a,int n) { DFT(a,n,inv(3)); int u=inv(n); for (int i=0;i<n;i++) a[i]=1ll*a[i]*u%P; } void mul(int *a,int *b,int n) { DFT(a,n,3),DFT(b,n,3); for (int i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%P; IDFT(a,n); } void Inv(int *a,int *b,int n) { if (n==1) {for (int i=0;i<t;i++) b[i]=0;b[0]=inv(a[0]);return;} Inv(a,b,n>>1); for (int i=0;i<n;i++) A[i]=a[i]; for (int i=n;i<(n<<1);i++) A[i]=0; n<<=1; DFT(A,n,3),DFT(b,n,3); for (int i=0;i<n;i++) b[i]=1ll*b[i]*(P+2-1ll*A[i]*b[i]%P)%P; IDFT(b,n); n>>=1; for (int i=n;i<(n<<1);i++) b[i]=0; } void trans(int *a,int *b,int n){for (int i=0;i<n;i++) b[i]=1ll*a[i+1]*(i+1)%P;} void dx(int *a,int *b,int n){b[0]=0;for (int i=1;i<n;i++) b[i]=1ll*a[i-1]*inv(i)%P;} void Ln(int *a,int t) { memset(x,0,sizeof(x)),memset(y,0,sizeof(y)); trans(a,x,t);Inv(a,y,t>>1);mul(x,y,t);dx(x,a,t); } void Exp(int *a,int *b,int n) { if (n==1){for (int i=0;i<t;i++) b[i]=0;b[0]=1;return;} Exp(a,b,n>>1); for (int i=0;i<(n>>1);i++) B[i]=b[i]; for (int i=(n>>1);i<n;i++) B[i]=0; Ln(B,n); for (int i=0;i<n;i++) B[i]=(P-B[i]+a[i])%P; B[0]=(B[0]+1)%P; for (int i=n;i<(n<<1);i++) B[i]=0; mul(b,B,n<<1); for (int i=n;i<(n<<1);i++) b[i]=0; } int main() { #ifndef ONLINE_JUDGE freopen("bag.in","r",stdin); freopen("bag.out","w",stdout); const char LL[]="%I64d\n"; #else const char LL[]="%lld\n"; #endif n=read(),m=read(); for (int i=1;i<=n;i++) a[read()]++;n=m; for (int i=1;i<=n;i++) for (int j=i;j<=n;j+=i) f[j]=(f[j]+1ll*a[i]*inv(j/i))%P; int t=1;while (t<=(n<<1)) t<<=1; Exp(f,g,t); for (int i=1;i<=m;i++) printf("%d\n",g[i]); return 0; }