洛谷 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) = sumlimits_{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 xⁿ ) 的话，可以这样算：

( F(x)=F_0(x)-frac{ G(F_0(x)) }{ G'(F_0(x)) } ) ( mod xⁿ ) ，其中 ( G(F_0(x)) equiv 0 (mod x^{ leftlceilfrac{n}{2} ight ceil } ) )

关于导数：

复合函数的导数：( 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 回去。相加就是埃氏筛的复杂度了。

转化就是这样：

( prodlimits_{i}frac{1}{1-x^{v_i}} => e^{sumlimits_{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})sumlimits_{j=1}^{infty}j*v_{i}*x^{j*v_{i}-1}dx )　　//变成从1开始了

　　　　　　　　　　(= int ( sumlimits_{j=1}^{infty}j*v_{i}*x^{j*v_{i}-1} - sumlimits_{j=1}^{infty}j*v_{i}*x^{(j+1)v_{i}-1} ) dx )

　　　　　　　　　　(= int sumlimits_{j=1}^{infty}v_{i}*x^{j*v_{i}-1} dx )

　　　　　　　　　　(= sumlimits_{j=1}^{infty}frac{1}{j}x^{j*v_{i}} )

注意相加的时候不要枚举每个物品，那样可能 n² ，应该枚举每种价值。

#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
",b[i]);
  return 0;
}