BZOJ 4913 [Sdoi2017] 遗忘的集合

题目链接

https://lydsy.com/JudgeOnline/problem.php?id=4913

题解

令$a_i=0/1$表示元素$i$是否在集合中，那么元素$i$的生成函数为
$(\frac{1}{1-x^i})^{a_i}$
现在已知了
$F(x)=\prod_{i=1}^{\infin}(\frac{1}{1-x^i})^{a_i}$
要求$a_i$。

对上式两边取对数
$\ln F(x)=\sum_{i=1}^{\infin}-a_i\ln(1-x^i)$
对右边泰勒展开
$\ln F(x)=\sum_{i=1}^{\infin}a_i\sum_{j=1}^{\infin}\frac{x^{ij}}{j}$
枚举$T=ij$
$\ln F(x)=\sum_{T=1}^{\infin}x^{T}\sum_{d|T}a_d\frac{d}{T}$
即
$T[x^T]\ln F(x)=\sum_{d|T}a_dd$
反演得到
$a_i=\frac{\sum_{d|i}\mu(\frac{i}{d})d[x^d]\ln F(x)}{i}$
因此只要多项式求$\ln$即可得到$a_i$的值。

注意常数优化，可以使用这篇文章提到的卡常技巧。

代码

#include <cmath>
#include <cstdio>
#include <algorithm>

#define getchar()(p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)

char buf[1<<21],*p1=buf,*p2=buf; 

int read()
{
  int x=0,f=1;
  char ch=getchar();
  while((ch<'0')||(ch>'9'))
    {
      if(ch=='-')
        {
          f=-f;
        }
      ch=getchar();
    }
  while((ch>='0')&&(ch<='9'))
    {
      x=x*10+ch-'0';
      ch=getchar();
    }
  return x*f;
}

int write(int x)
{
  if(x>9)
    {
      write(x/10);
    }
  putchar(x%10+'0');
  return 0;
}

const int maxn=524288;
const double pi=acos(-1);

struct complex
{
  double r,i;

  complex(double _r=0,double _i=0):r(_r),i(_i){}

  complex operator +(const complex &oth) const
  {
    return complex(r+oth.r,i+oth.i);
  }

  complex operator -(const complex &oth) const
  {
    return complex(r-oth.r,i-oth.i);
  }

  complex operator *(const complex &oth) const
  {
    return complex(r*oth.r-i*oth.i,r*oth.i+i*oth.r);
  }
};

inline complex conj(complex x)
{
  return complex(x.r,-x.i);
}

const complex hreal=complex(0.5,0);
const complex himag=complex(0,-0.5);
const complex gimag=complex(0,1);

complex w[maxn+10];

inline int fft(complex *s,int len)
{
  for(int i=0,j=0; i<len; ++i)
    {
      if(j<i)
        {
          std::swap(s[j],s[i]);
        }
      int l=len>>1;
      while((j^=l)<l)
        {
          l>>=1;
        }
    }
  for(int i=1,d=maxn; i<len; i<<=1,d>>=1)
    {
      for(int j=0; j<len; j+=(i<<1))
        {
          for(int k=j; k<j+i; ++k)
            {
              complex x=s[k],y=w[d*(k-j)]*s[k+i];
              s[k]=x+y;
              s[k+i]=x-y;
            }
        }
    }
  return 0;
}

int mod;
complex p[maxn+10],q[maxn+10],da[maxn+10],db[maxn+10],dc[maxn+10];

int multi(int *a,int la,int *b,int lb,int *res,int mx)
{
  int l=0,n=1;
  while(n<la+lb-1)
    {
      ++l;
      n<<=1;
    }
  for(int i=0; i<la; ++i)
    {
      p[i]=complex(a[i]>>15,a[i]&32767);
    }
  for(int i=0; i<lb; ++i)
    {
      q[i]=complex(b[i]>>15,b[i]&32767);
    }
  for(int i=la; i<n; ++i)
    {
      p[i]=complex(0,0);
    }
  for(int i=lb; i<n; ++i)
    {
      q[i]=complex(0,0);
    }
  fft(p,n);
  fft(q,n);
  for(int i=0; i<n; ++i)
    {
      int j=(n-i)&(n-1);
      complex tka=hreal*(p[i]+conj(p[j])),
        tra=himag*(p[i]-conj(p[j])),
        tkb=hreal*(q[i]+conj(q[j])),
        trb=himag*(q[i]-conj(q[j]));
      da[j]=tka*tkb;
      db[j]=tra*trb;
      dc[j]=tka*trb+tra*tkb;
    }
  for(int i=0; i<n; ++i)
    {
      p[i]=da[i]+gimag*db[i];
    }
  fft(p,n);
  fft(dc,n);
  for(int i=0; i<mx; ++i)
    {
      long long ac=(long long)(p[i].r/n+0.5),
        bd=(long long)(p[i].i/n+0.5),
        ad=(long long)(dc[i].r/n+0.5);
      res[i]=(bd+((ad%mod)<<15)+((ac%mod)<<30))%mod;
    }
  return 0;
}

inline int quickpow(int a,int b)
{
  int res=1;
  while(b)
    {
      if(b&1)
        {
          res=1ll*res*a%mod;
        }
      a=1ll*a*a%mod;
      b>>=1;
    }
  return res;
}

int tmp[maxn+10];

inline int getinv(int *s,int len,int *res)
{
  if(len==1)
    {
      res[0]=quickpow(s[0],mod-2);
      return 0;
    }
  int k=(len+1)/2;
  getinv(s,k,res);
  for(int i=0; i<len; ++i)
    {
      tmp[i]=s[i];
    }
  multi(res,k,tmp,len,tmp,len);
  multi(res,k,tmp+k,len-k,res+k,len-k);
  for(int i=k; i<len; ++i)
    {
      if(res[i])
        {
          res[i]=mod-res[i];
        }
    }
  return 0;
}

inline int der(int *s,int len,int *res)
{
  for(int i=0; i<len-1; ++i)
    {
      res[i]=1ll*s[i+1]*(i+1)%mod;
    }
  res[len-1]=0;
  return 0;
}

int inv[maxn+10];

inline int inter(int *s,int len,int *res)
{
  for(int i=len; i; --i)
    {
      res[i]=1ll*s[i-1]*inv[i]%mod;
    }
  res[0]=0;
  return 0;
}

int ln[maxn+10];

int getln(int *s,int len,int *res)
{
  der(s,len,ln);
  getinv(s,len,res);
  multi(ln,len-1,res,len,res,len-1);
  inter(res,len-1,res);
  return 0;
}

int n,f[maxn+10],ip[maxn+10],prime[maxn+10],cnt,mu[maxn+10],g[maxn+10],r[maxn+10];

int getprime()
{
  ip[1]=mu[1]=1;
  for(int i=2; i<=maxn; ++i)
    {
      if(!ip[i])
        {
          prime[++cnt]=i;
          mu[i]=-1;
        }
      for(int j=1; (j<=cnt)&&(i*prime[j]<=maxn); ++j)
        {
          int x=i*prime[j];
          ip[x]=1;
          if(i%prime[j]==0)
            {
              mu[x]=0;
              break;
            }
          mu[x]=-mu[i];
        }
    }
  double k=pi/maxn;
  for(int i=0; i<maxn; ++i)
    {
      w[i]=complex(cos(i*k),sin(i*k));
    }
  inv[0]=inv[1]=1;
  for(int i=2; i<=maxn; ++i)
    {
      inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    }
  return 0;
}

int ans[maxn+10];

int main()
{
  n=read();
  mod=read();
  getprime();
  r[0]=1;
  for(int i=1; i<=n; ++i)
    {
      r[i]=read();
    }
  getln(r,n+1,f);
  for(int i=1; i<=n; ++i)
    {
      int x=1ll*i*f[i]%mod;
      for(int j=i; j<=n; j+=i)
        {
          g[j]=g[j]+mu[j/i]*x;
          if(g[j]<0)
            {
              g[j]+=mod;
            }
          if(g[j]>=mod)
            {
              g[j]-=mod;
            }
        }
    }
  int tot=0;
  for(int i=1; i<=n; ++i)
    {
      if(g[i])
        {
          ans[++tot]=i;
        }
    }
  write(tot);
  putchar('\n');
  for(int i=1; i<tot; ++i)
    {
      write(ans[i]);
      putchar(' ');
    }
  write(ans[tot]);
  putchar('\n');
  return 0;
}