• 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
    ";
    #else
    	const char LL[]="%lld
    ";
    #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
    ",g[i]);
    	return 0;
    }
    

      

  • 相关阅读:
    终于看到费德勒在法网如愿!
    o(∩_∩)o...,今天去博客园了!
    条款4:使用Conditional特性代替#if条件编译
    MSDTC无法启动的解决方法
    2009 很有意义的一天
    从现在开始,争取记录每天所学到的、所感受到的、所遇见到的点点滴滴!
    了解MOSS2007 内容类型ID(Content Type IDs)命名规则
    CreateSpecificCulture('zhcn')和new CultureInfo('zhcn')的区别
    金华大显数码科技有限公司诚聘
    使用SQL Server中按位于来表示组合状态
  • 原文地址:https://www.cnblogs.com/Gloid/p/10393137.html
Copyright © 2020-2023  润新知