[luogu4389]付公主的背包

前言

一道巧妙的推式子题

题目相关

链接

题目大意

给出$n$个商品，第$i$个商品的体积为viv_ivi​，并且有无限个
现在给出一个$m$，对背包大小$s\in [1,m]$求填满方案数

数据范围

$n\le 100000,m\le 100000$

做法

算法一

考虑生成函数
对于第$i$个物品我们构造生成函数fi(x)=∑j=0∞xvijf_i(x)=\sum_{j=0}^{\infty}x^{v_ij}fi​(x)=j=0∑∞​xvi​j
所有生成函数相乘∏i=1nfi(x)\prod_{i=1}^nf_i(x)i=1∏n​fi​(x)
显然答案的$[1,m]$次系数就是答案了
我们发现整个算法的复杂度是$\mathcal{O}(nmlogm)$，显然不能通过

算法二

我们发现生成函数其实可以看成等比数列，利用求和公式
我们推一下
fi(x)=∑j=0∞xvijf_i(x)=\sum_{j=0}^{\infty}x^{v_ij}fi​(x)=j=0∑∞​xvi​j
xvifi(x)=∑j=1∞xvijx^{v_i}f_i(x)=\sum_{j=1}^{\infty}x^{v_ij}xvi​fi​(x)=j=1∑∞​xvi​j
相减得fi(x)−xvifi(x)=1f_i(x)-x^{v_i}f_i(x)=1fi​(x)−xvi​fi​(x)=1
(1−xvi)fi(x)=1(1-x^{v_i})f_i(x)=1(1−xvi​)fi​(x)=1
即fi(x)=11−xvif_i(x)=\frac1{1-x^{v_i}}fi​(x)=1−xvi​1​

算法三

现在我们要快速求∏i=1n1(1−xvi)\prod_{i=1}^n\frac1{(1-x^{v_i})}i=1∏n​(1−xvi​)1​
一个很好的思路是$ln$后进行加法操作以完成原式的乘，再$exp$即可
那么现在要求ln(1(1−xvi))ln(\frac1{(1-x^{v_i})})ln((1−xvi​)1​)

让我们来推一下：
我们发现，直接做并不好做
考虑回到原式求导后积分
我们已知A(x)=11−xviA(x)=\frac1{1-x^{v_i}}A(x)=1−xvi​1​，设$B(x)=ln(A(x))$
我们要求B′(x)B'(x)B′(x)的值的话先把A′(x)A'(x)A′(x)求出
我们回最上面的式子进行求导可知
A′(x)=(∑j=0∞xvij)′=∑j=0∞(xvij)′=∑j=0∞vijxvij−1\begin{aligned} A'(x)&=(\sum_{j=0}^{\infty}x^{v_ij})'\\ &=\sum_{j=0}^{\infty}(x^{v_ij})'\\ &=\sum_{j=0}^{\infty}v_ijx^{v_ij-1} \end{aligned}A′(x)​=(j=0∑∞​xvi​j)′=j=0∑∞​(xvi​j)′=j=0∑∞​vi​jxvi​j−1​
（之后写的时候省略$i$）
根据复合函数的求导法则可知
B′(x)=A′(x)A(x)=∑j=0∞vijxvj−111−xv=(1−xv)∑j=0∞vjxvj−1=(∑j=0∞vjxvj−1)−(∑j=0∞vjxv(j+1)−1)=(∑j=0∞vjxvj−1)−(∑j=1∞v(j−1)xvj−1)=(∑j=1∞vjxvj−1)−(∑j=1∞v(j−1)xvj−1)=∑j=1∞v(j−(j−1))xvj−1=∑j=1∞vxvj−1\begin{aligned} B'(x)&=\frac{A'(x)}{A(x)}\\ &=\frac{\sum_{j=0}^{\infty}v_ijx^{vj-1}}{\frac1{1-x^v}}\\ &=(1-x^v)\sum_{j=0}^{\infty}vjx^{vj-1}\\ &=\left(\sum_{j=0}^{\infty}vjx^{vj-1}\right)-\left(\sum_{j=0}^{\infty}vjx^{v(j+1)-1}\right)\\ &=\left(\sum_{j=0}^{\infty}vjx^{vj-1}\right)-\left(\sum_{j=1}^{\infty}v(j-1)x^{vj-1}\right)\\ &=\left(\sum_{j=1}^{\infty}vjx^{vj-1}\right)-\left(\sum_{j=1}^{\infty}v(j-1)x^{vj-1}\right)\\ &=\sum_{j=1}^{\infty}v(j-(j-1))x^{vj-1}\\ &=\sum_{j=1}^{\infty}vx^{vj-1}\\ \end{aligned}B′(x)​=A(x)A′(x)​=1−xv1​∑j=0∞​vi​jxvj−1​=(1−xv)j=0∑∞​vjxvj−1=(j=0∑∞​vjxvj−1)−(j=0∑∞​vjxv(j+1)−1)=(j=0∑∞​vjxvj−1)−(j=1∑∞​v(j−1)xvj−1)=(j=1∑∞​vjxvj−1)−(j=1∑∞​v(j−1)xvj−1)=j=1∑∞​v(j−(j−1))xvj−1=j=1∑∞​vxvj−1​
然后我们重新积分
∫B′(x)=∫(∑j=1∞vxvj−1)=∑j=1∞∫vxvj−1=∑j=1∞vxvjvj=∑j=1∞1jxvj\begin{aligned} \int B'(x)&=\int \left(\sum_{j=1}^{\infty}vx^{vj-1}\right)\\ &=\sum_{j=1}^{\infty}\int vx^{vj-1}\\ &=\sum_{j=1}^{\infty}v\frac{x^{vj}}{vj}\\ &=\sum_{j=1}^{\infty}\frac1jx^{vj}\\ \end{aligned}∫B′(x)​=∫(j=1∑∞​vxvj−1)=j=1∑∞​∫vxvj−1=j=1∑∞​vvjxvj​=j=1∑∞​j1​xvj​

我们现在需要快速求出∑i=1n∑j=1∞1jxvij\sum_{i=1}^n\sum_{j=1}^{\infty}\frac1jx^{v_ij}i=1∑n​j=1∑∞​j1​xvi​j
我们发现，后者的次数是vijv_ijvi​j的形式的，所以我们将所有相同的viv_ivi​合并，那么我们就可以通过调和级数$\mathcal{O}(mlogm)$的时间复杂度获得整个数组
然后$exp$即可，复杂度$\mathcal{O}(nlogn)$

代码

#include
#include
#include
#include
#include
namespace fast_IO
{const int IN_LEN=10000000,OUT_LEN=10000000;char ibuf[IN_LEN],obuf[OUT_LEN],*ih=ibuf+IN_LEN,*oh=obuf,*lastin=ibuf+IN_LEN,*lastout=obuf+OUT_LEN-1;inline char getchar_(){return (ih==lastin)&&(lastin=(ih=ibuf)+fread(ibuf,1,IN_LEN,stdin),ih==lastin)?EOF:*ih++;}inline void putchar_(const char x){if(oh==lastout)fwrite(obuf,1,oh-obuf,stdout),oh=obuf;*oh++=x;}inline void flush(){fwrite(obuf,1,oh-obuf,stdout);}
}
using namespace fast_IO;
#define getchar() getchar_()
#define putchar(x) putchar_((x))
typedef long long ll;
#define rg register
template <typename T> inline T max(const T a,const T b){return a>b?a:b;}
template <typename T> inline T min(const T a,const T b){return a<b?a:b;}
template <typename T> inline T mind(T&a,const T b){a=a<b?a:b;}
template <typename T> inline T maxd(T&a,const T b){a=a>b?a:b;}
template <typename T> inline T abs(const T a){return a>0?a:-a;}
template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;}
template <typename T> inline void swap(T*a,T*b){T c=a;a=b;b=c;}
template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);}
template <typename T> inline T square(const T x){return x*x;};
template <typename T> inline void read(T&x)
{char cu=getchar();x=0;bool fla=0;while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}while(isdigit(cu))x=x*10+cu-'0',cu=getchar();if(fla)x=-x;  
}
template <typename T> void printe(const T x)
{if(x>=10)printe(x/10);putchar(x%10+'0');
}
template <typename T> inline void print(const T x)
{if(x<0)putchar('-'),printe(-x);else printe(x);
}
const int maxn=2097152,mod=998244353;
inline int Md(const int x){return x>=mod?x-mod:x;}
template<typename T>
inline int pow(int x,T y)
{rg int res=1;x%=mod;for(;y;y>>=1,x=(ll)x*x%mod)if(y&1)res=(ll)res*x%mod;return res;
}
namespace Poly///namespace of Poly
{
int W_[maxn],ha[maxn],hb[maxn],Inv[maxn];
inline void init(const int x)
{rg int tim=0,lenth=1;while(lenth<x)lenth<<=1,tim++;for(rg int i=1;i<lenth;i<<=1){const int WW=pow(3,(mod-1)/(i*2));W_[i]=1;for(rg int j=i+1,k=i<<1;j<k;j++)W_[j]=(ll)W_[j-1]*WW%mod;}Inv[0]=Inv[1]=1;for(rg int i=2;i<x;i++)Inv[i]=(ll)(mod-mod/i)*Inv[mod%i]%mod;
}
int L;
inline void DFT(int*A)//prepare:init L 
{for(rg int i=0,j=0;i<L;i++){if(i>j)swap(A[i],A[j]);for(rg int k=L>>1;(j^=k)<k;k>>=1);}for(rg int i=1;i<L;i<<=1)for(rg int j=0,J=i<<1;j<L;j+=J)for(rg int k=0;k<i;k++){const int x=A[j+k],y=(ll)A[j+k+i]*W_[i+k]%mod;A[j+k]=Md(x+y),A[j+k+i]=Md(mod+x-y);}
}
void IDFT(int*A)
{for(rg int i=1;i<L-i;i++)swap(A[i],A[L-i]);DFT(A);
}
inline int Quadratic_residue(const int a)
{if(a==0)return 0;int b=(rand()<<14^rand())%mod;while(pow(b,(mod-1)>>1)!=mod-1)b=(rand()<<14^rand())%mod;int s=mod-1,t=0,x,inv=pow(a,mod-2),f=1;while(!(s&1))s>>=1,t++,f<<=1;t--,x=pow(a,(s+1)>>1),f>>=1;while(t){f>>=1;if(pow((int)((ll)inv*x%mod*x%mod),f)!=1)x=(ll)x*pow(b,s)%mod;t--,s<<=1;}return min(x,mod-x);
}
struct poly
{std::vector<int>A;poly(){A.resize(0);}poly(const int x){A.resize(1),A[0]=x;}inline int&operator[](const int x){return A[x];}inline int operator[](const int x)const{return A[x];}inline void clear(){A.clear();}inline unsigned int size()const{return A.size();}inline void resize(const unsigned int x){A.resize(x);}void RE(const int x){A.resize(x);for(rg int i=0;i<x;i++)A[i]=0; }void readin(const int MAX){A.resize(MAX);for(rg int i=0;i<MAX;i++)read(A[i]);}void putout()const{for(rg unsigned int i=0;i<A.size();i++)print(A[i]),putchar(' ');}inline poly operator +(const poly b)const{poly RES;RES.resize(max(size(),b.size()));for(rg unsigned int i=0;i<RES.size();i++)RES[i]=Md((i<size()?A[i]:0)+(i<b.size()?b[i]:0));return RES;}inline poly operator -(const poly b)const{poly RES;RES.resize(max(size(),b.size()));for(rg unsigned int i=0;i<RES.size();i++)RES[i]=Md((i<size()?A[i]:0)+mod-(i<b.size()?b[i]:0));return RES;}inline poly operator *(const int b)const{poly RES=*this;for(rg unsigned int i=0;i<RES.size();i++)RES[i]=(ll)RES[i]*b%mod;return RES;}inline poly operator /(const int b)const{poly RES=(*this)*pow(b,mod-2);return RES;}inline poly operator *(const poly b)const{const int RES=A.size()+b.size()-1;L=1;while(L<RES)L<<=1;poly c;c.A.resize(RES);memset(ha,0,L<<2);memset(hb,0,L<<2);for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];for(rg unsigned int i=0;i<b.A.size();i++)hb[i]=b.A[i];DFT(ha),DFT(hb);for(rg int i=0;i<L;i++)ha[i]=(ll)ha[i]*hb[i]%mod;IDFT(ha);const int inv=pow(L,mod-2);for(rg int i=0;i<RES;i++)c.A[i]=(ll)ha[i]*inv%mod;return c;}inline poly inv()const{poly C;if(A.size()==1){C=*this;C[0]=pow(C[0],mod-2);return C;}C.resize((A.size()+1)>>1);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];C=C.inv();L=1;while(L<(int)size()*2-1)L<<=1;for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];for(rg unsigned int i=0;i<C.size();i++)hb[i]=C[i];memset(ha+A.size(),0,(L-A.size())<<2);memset(hb+C.size(),0,(L-C.size())<<2);DFT(ha),DFT(hb);for(rg int i=0;i<L;i++)ha[i]=(ll)(2-(ll)hb[i]*ha[i]%mod+mod)*hb[i]%mod;IDFT(ha);const int inv=pow(L,mod-2);C.resize(size());for(rg unsigned int i=0;i<size();i++)C[i]=(ll)ha[i]*inv%mod;return C;}
/*    inline poly inv()const{poly C;if(A.size()==1){C=*this;C[0]=pow(C[0],mod-2);return C;}C.resize((A.size()+1)>>1);for(rg unsigned int i=0;i//大常数版本 inline void Reverse(const int n){A.resize(n);for(rg int i=0,j=n-1;i<j;i++,j--)swap(A[i],A[j]);}inline poly operator /(const poly B)const{if(size()<B.size())return 0;poly a=*this,b=B;a.Reverse(size()),b.Reverse(B.size());b.resize(size()-B.size()+1);b=b.inv();b=b*a;b.Reverse(size()-B.size()+1);return b;}inline poly operator %(const poly B)const{poly C=(*this)-(*this)/B*B;C.resize(B.size()-1);return C;}inline poly diff()const{poly C;C.resize(size()-1);for(rg unsigned int i=1;i<size();i++)C[i-1]=(ll)A[i]*i%mod;return C;}inline poly inte()const{poly C;C.resize(size()+1);for(rg unsigned int i=0;i<size();i++)C[i+1]=(ll)A[i]*Inv[i+1]%mod;C[0]=0;return C;}inline poly ln()const{poly C=(diff()*inv()).inte();C.resize(size());return C;}inline poly exp()const{poly C;if(size()==1){C=*this;C[0]=1;return C;}C.resize((size()+1)>>1);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];C=C.exp();C.resize(size());poly D=(poly)1-C.ln()+*this;D=D*C;D.resize(size());return D;}inline poly sqrt()const{poly C;if(size()==1){C=*this;C[0]=Quadratic_residue(C[0]);return C;}C.resize((size()+1)>>1);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];C=C.sqrt();C.resize(size());C=(C+*this*C.inv())*(int)499122177;C.resize(size());return C;}inline poly operator >>(const unsigned int x)const{poly C;if(size()<x){C.resize(0);return C;}C.resize(size()-x);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i+x];return C;}inline poly operator <<(const unsigned int x)const{poly C;C.RE(size()+x);for(rg unsigned int i=0;i<size();i++)C[i+x]=A[i];return C;}inline poly Pow(const unsigned int x)const{for(rg unsigned int i=0;i<size();i++)if(A[i]){poly C=((((*this/A[i])>>i).ln()*x).exp()*pow(A[i],x))<<(min(i*x,size()));C.resize(size());return C;}return *this;}inline void cheng(const poly&B){for(rg unsigned int i=0;i<size();i++)A[i]=(ll)A[i]*B[i]%mod; }inline void jia(const poly&B){for(rg unsigned int i=0;i<size();i++)A[i]=Md(A[i]+B[i]); }inline void dft(){memset(ha,0,L<<2);for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];DFT(ha);resize(L);for(rg int i=0;i<L;i++)A[i]=ha[i];}inline void idft(){memset(ha,0,L<<2);for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];IDFT(ha);const int inv=pow(L,mod-2);for(rg int i=0;i<L;i++)A[i]=(ll)ha[i]*inv%mod;while(size()&&!A[size()-1])A.pop_back();}
};
void fz(const int root,const int l,const int r,std::vector<int>&v,std::vector<poly>&A)
{if(l==r){A[root].resize(2);A[root][0]=(mod-v[l])%mod;A[root][1]=1;return;}const int mid=(l+r)>>1;fz(root<<1,l,mid,v,A),fz(root<<1|1,mid+1,r,v,A);A[root]=A[root<<1]*A[root<<1|1];
}
void calc(const int root,const int l,const int r,std::vector<int>&v,std::vector<poly>&A,std::vector<poly>&B)
{if(l==r){v[l]=B[root][0];return;}const int mid=(l+r)>>1;B[root<<1]=B[root]%A[root<<1];B[root<<1|1]=B[root]%A[root<<1|1];calc(root<<1,l,mid,v,A,B),calc(root<<1|1,mid+1,r,v,A,B);
}
void multi_point_evaluation(const poly a,std::vector<int>&v)
{std::vector<poly>A,B;A.resize(maxn),B.resize(maxn);fz(1,0,v.size()-1,v,A);B[1]=a%A[1];calc(1,0,v.size()-1,v,A,B);
}
void fz2(const int root,const int l,const int r,std::vector<int>&y,std::vector<poly>&A,std::vector<poly>&B)
{if(l==r){B[root].resize(1),B[root][0]=y[l];return;}const int mid=(l+r)>>1;fz2(root<<1,l,mid,y,A,B),fz2(root<<1|1,mid+1,r,y,A,B);B[root]=B[root<<1]*A[root<<1|1]+B[root<<1|1]*A[root<<1];
}
poly interpolation(std::vector<int>&x,std::vector<int>&y)
{std::vector<poly>A,B;A.resize(maxn),B.resize(maxn);fz(1,0,x.size()-1,x,A);multi_point_evaluation(A[1].diff(),x);for(rg unsigned int i=0;i<x.size();i++)y[i]=(ll)y[i]*pow(x[i],mod-2)%mod;fz2(1,0,x.size()-1,y,A,B);return B[1];
}
}///namespace of Poly
int n,m;Poly::poly a;
int main()
{Poly::init(maxn);///namespace of Polyread(n),read(m);a.resize(m+1);for(rg int i=1;i<=n;i++){int v;read(v);a[v]++;}for(rg int i=m;i>=1;i--)for(rg int j=2;i*j<=m;j++)a[i*j]=Md((ll)a[i]*(Poly::Inv[j])%mod+a[i*j]);a=a.exp();for(rg int i=1;i<=m;i++)print(a[i]),putchar('\n');return flush(),0;
}

总结

还有一道完全加强版的题目loj556
之后我会更博客

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！