题目: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=;bool fx=;char ch=getchar();
while(ch>''||ch<''){if(ch=='-')fx=;ch=getchar();}
while(ch>=''&&ch<='')ret=ret*+ch-'',ch=getchar();
return fx?ret:-ret;
}
int Mx(int a,int b){return a>b?a:b;}
const int N=(<<)+,mod=,gen=;
int upt(int x){while(x>=mod)x-=mod;while(x<)x+=mod;return x;}
int pw(int x,int k)
{int ret=;while(k){if(k&)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=;}return ret;} int n,m,cnt[N],a[N],b[N];
namespace poly{
int len,r[N],inv[N],Wn[N][],tp[N],A[N],B[N],tp2[N];
void init(int len)
{
inv[]=;
for(int i=;i<=len;i++)
inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod;
for(int R=;R<=len;R<<=)
Wn[R][]=pw( ,(mod-)/R ),
Wn[R][]=pw( ,(mod-)-(mod-)/R );
}
void ntt_pre(int len)
{
for(int i=,j=len>>;i<len;i++)
r[i]=(r[i>>]>>)+((i&)?j:);
}
void ntt(int *a,bool fx)
{
for(int i=;i<len;i++)
if(i<r[i])swap(a[i],a[r[i]]);
for(int R=;R<=len;R<<=)
{
int wn=Wn[R][fx];
for(int i=,m=R>>;i<len;i+=R)
for(int j=,w=;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=;i<len;i++)a[i]=(ll)a[i]*t%mod;
}
void get_dao(int n,int *a,int *b)
{
for(int i=;i<n;i++)b[i-]=(ll)a[i]*i%mod;
b[n-]=;
}
void get_jf(int n,int *a,int *b)
{
for(int i=;i<n;i++)b[i]=(ll)a[i-]*inv[i]%mod;
b[]=;
}
void get_inv(int n,int *a,int *b)
{
b[]=pw(a[],mod-);
for(int tn=,l=;tn<n;tn=l,l<<=)
{
for(int i=tn;i<l;i++)b[i]=;
for(int i=;i<l;i++)tp2[i]=a[i],tp2[i+l]=;
len=l<<; ntt_pre(len);
ntt(tp2,); ntt(b,);
for(int i=;i<len;i++)
b[i]=(ll)b[i]*(-(ll)tp2[i]*b[i]%mod+mod)%mod;
ntt(b,);
}
}
void get_ln(int n,int *a,int *b)
{
for(int i=;i<n;i++)B[i]=;///
get_dao(n,a,A);get_inv(n,a,B);
len=n<<; ntt_pre(len);
ntt(A,); ntt(B,);
for(int i=;i<len;i++)A[i]=(ll)A[i]*B[i]%mod;
ntt(A,); get_jf(n,A,b);
}
void get_exp(int n,int *a,int *b)
{
b[]=;
for(int tn=,l=;tn<n;tn=l,l<<=)//b[0~tn-1]->b[0~l-1]
{
for(int i=tn;i<l;i++)b[i]=; get_ln(l,b,tp);
for(int i=;i<l;i++)
tp[i]=upt(-tp[i]+a[i]); tp[]=upt(tp[]+);
len=l<<; ntt_pre(len);
ntt(b,); ntt(tp,);
for(int i=;i<len;i++)b[i]=(ll)tp[i]*b[i]%mod;
ntt(b,);
}
}
}
int main()
{
n=rdn();m=rdn(); int l;for(l=;l<=m;l<<=);//m not n!
poly::init(l<<);//l<<1
for(int i=,d;i<=n;i++)d=rdn(),cnt[d]++;
for(int i=;i<=m;i++)
{
if(!cnt[i])continue;
for(int j=;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=;i<=m;i++)printf("%d\n",b[i]);
return ;
}