洛谷 4389 付公主的背包——多项式求ln、exp
题目:https://www.luogu.org/problemnew/show/P4389
关于泰勒展开:
https://blog.csdn.net/SoHardToNamed/article/details/80550935
https://www.cnblogs.com/guo-xiang/p/6662881.html
大概就是:\( f(x) = \sum\limits_{i=0}^{n}\frac{ f^{(i)}(x_0) }{i!}(x-x_0)^i +R_n\)
麦克劳林展开就是上面的 \( x_0 \) 取 0 的时候。
关于牛顿迭代:
https://www.matongxue.com/madocs/205.html
//blog.miskcoo.com/2015/06/polynomial-with-newton-method
大概就是如果有 \( F(G(x)) \equiv 0 \) ( mod xn ) 的话,可以这样算:
\( F(x)=F_0(x)-\frac{ G(F_0(x)) }{ G'(F_0(x)) } \) ( mod xn ) ,其中 \( G(F_0(x)) \equiv 0 (mod x^{ \left\lceil\frac{n}{2}\right\rceil } ) \)
关于导数:
复合函数的导数:( F( G( x ) ) )' = F'( G( x ) ) * G'( x )
分数的导数:\( \frac{U}{V} = \frac{U'V-UV'}{V^2} \)
所以求 exp 就是这样的(参见上面粘的 miskcoo 的博客):
\( f(x) = e^{A(x)} \) 即 \( ln(f(x))-A(x) = 0 \)
令 \( G(f(x)) = ln(f(x))-A(x) \) ,则 \( G'(f(x)) = \frac{1}{f(x)} \)
\( f(x) = f_0(x)-\frac{G(f_0(x))}{G'(f_0(x))} \)
\(= f_0(x)-\frac{ ln(f_0(x))-A(x) }{ \frac{1}{f_0(x)} } \)
\(= f_0(x)(1-ln(f_0(x))+A(x)) \)
关于本题,就是写出每种物品的生成函数,需要把它们都乘起来,但可以转化成 ln 相加再 exp 回去。相加就是埃氏筛的复杂度了。
转化就是这样:
\( \prod\limits_{i}\frac{1}{1-x^{v_i}} => e^{\sum\limits_{i}ln( \frac{1}{1-x^{v_i}} )} \)
\( ln(\frac{1}{1-x^{v_i}}) = \int ( ln(\frac{1}{1-x^{v_i}}) )' dx \)
这里有一个套路:\( ( ln(f(x)) )' = \frac{f'(x)}{f(x)} \)
现在想让一会儿积分回去的式子是一个比较好求的式子,一般是一个 \( \sum \) 样子的式子。
所以把分子写成 \( \sum \) 的样子,然后把分母代入每一项里,就能得到一个新的 \( \sum \) 的样子,一般就可以了。
\(= \int (1-x^{v_i})\sum\limits_{j=1}^{\infty}j*v_{i}*x^{j*v_{i}-1}dx \) //变成从1开始了
\(= \int ( \sum\limits_{j=1}^{\infty}j*v_{i}*x^{j*v_{i}-1} - \sum\limits_{j=1}^{\infty}j*v_{i}*x^{(j+1)v_{i}-1} ) dx \)
\(= \int \sum\limits_{j=1}^{\infty}v_{i}*x^{j*v_{i}-1} dx \)
\(= \sum\limits_{j=1}^{\infty}\frac{1}{j}x^{j*v_{i}} \)
注意相加的时候不要枚举每个物品,那样可能 n2 ,应该枚举每种价值。
#include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; int rdn() { int ret=0;bool fx=1;char ch=getchar(); while(ch>'9'||ch<'0'){if(ch=='-')fx=0;ch=getchar();} while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar(); return fx?ret:-ret; } int Mx(int a,int b){return a>b?a:b;} const int N=(1<<18)+5,mod=998244353,gen=3; int upt(int x){while(x>=mod)x-=mod;while(x<0)x+=mod;return x;} int pw(int x,int k) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;} int n,m,cnt[N],a[N],b[N]; namespace poly{ int len,r[N],inv[N],Wn[N][2],tp[N],A[N],B[N],tp2[N]; void init(int len) { inv[1]=1; for(int i=2;i<=len;i++) inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod; for(int R=2;R<=len;R<<=1) Wn[R][0]=pw( 3,(mod-1)/R ), Wn[R][1]=pw( 3,(mod-1)-(mod-1)/R ); } void ntt_pre(int len) { for(int i=0,j=len>>1;i<len;i++) r[i]=(r[i>>1]>>1)+((i&1)?j:0); } void ntt(int *a,bool fx) { for(int i=0;i<len;i++) if(i<r[i])swap(a[i],a[r[i]]); for(int R=2;R<=len;R<<=1) { int wn=Wn[R][fx]; for(int i=0,m=R>>1;i<len;i+=R) for(int j=0,w=1;j<m;j++,w=(ll)w*wn%mod) { int x=a[i+j], y=(ll)w*a[i+m+j]%mod; a[i+j]=upt(x+y); a[i+m+j]=upt(x-y); } } if(!fx)return; int t=inv[len]; for(int i=0;i<len;i++)a[i]=(ll)a[i]*t%mod; } void get_dao(int n,int *a,int *b) { for(int i=1;i<n;i++)b[i-1]=(ll)a[i]*i%mod; b[n-1]=0; } void get_jf(int n,int *a,int *b) { for(int i=1;i<n;i++)b[i]=(ll)a[i-1]*inv[i]%mod; b[0]=0; } void get_inv(int n,int *a,int *b) { b[0]=pw(a[0],mod-2); for(int tn=1,l=2;tn<n;tn=l,l<<=1) { for(int i=tn;i<l;i++)b[i]=0; for(int i=0;i<l;i++)tp2[i]=a[i],tp2[i+l]=0; len=l<<1; ntt_pre(len); ntt(tp2,0); ntt(b,0); for(int i=0;i<len;i++) b[i]=(ll)b[i]*(2-(ll)tp2[i]*b[i]%mod+mod)%mod; ntt(b,1); } } void get_ln(int n,int *a,int *b) { for(int i=0;i<n;i++)B[i]=0;/// get_dao(n,a,A);get_inv(n,a,B); len=n<<1; ntt_pre(len); ntt(A,0); ntt(B,0); for(int i=0;i<len;i++)A[i]=(ll)A[i]*B[i]%mod; ntt(A,1); get_jf(n,A,b); } void get_exp(int n,int *a,int *b) { b[0]=1; for(int tn=1,l=2;tn<n;tn=l,l<<=1)//b[0~tn-1]->b[0~l-1] { for(int i=tn;i<l;i++)b[i]=0; get_ln(l,b,tp); for(int i=0;i<l;i++) tp[i]=upt(-tp[i]+a[i]); tp[0]=upt(tp[0]+1); len=l<<1; ntt_pre(len); ntt(b,0); ntt(tp,0); for(int i=0;i<len;i++)b[i]=(ll)tp[i]*b[i]%mod; ntt(b,1); } } } int main() { n=rdn();m=rdn(); int l;for(l=1;l<=m;l<<=1);//m not n! poly::init(l<<1);//l<<1 for(int i=1,d;i<=n;i++)d=rdn(),cnt[d]++; for(int i=1;i<=m;i++) { if(!cnt[i])continue; for(int j=1;j*i<=m;j++) a[j*i]=(a[j*i]+(ll)cnt[i]*poly::inv[j])%mod; } poly::get_exp(l,a,b); for(int i=1;i<=m;i++)printf("%d\n",b[i]); return 0; }