bzoj 3456 城市规划 —— 分治FFT / 多项式求逆 / 指数型生成函数(多项式求ln)

题目：https://www.lydsy.com/JudgeOnline/problem.php?id=3456

首先考虑DP做法，正难则反，考虑所有情况减去不连通的情况；

而不连通的情况就是那个经典做法：选定一个划分点，枚举包含它的连通块，连通块以外的部分随便连(但不和连通块连通)，合起来就是不连通的方案数；

设 ( f[i] ) 表示一共 ( i ) 个点时的连通方案数，( g[i] ) 表示 ( i ) 个点随便连的方案数，即 ( g[i] = 2^{C_{i}^{2}} )，则：

( f[i] = 2^{C_{i}^{2}} - sumlimits_{j=1}^{i-1} C_{i-1}^{j-1} * f_{j} * g_{i-j} )

只要令 ( g[0] = 0 )，( j ) 就可以枚举到 ( i )；

把 ( C_{i-1}^{j-1} ) 拆开，对应分配到各处，得到：

( f[i] = 2^{C_{i}^{2}} - (i-1)!*sumlimits_{j=1}^{i}frac{f_{j}}{(j-1)!} * frac{g_{i-j}}{(i-j)!} )

所以把 ( g[i] ) 的定义改成 ( g[i] = frac{2^{C_{i}^{2}}}{i!} )

于是 ( f[i] = 2^{C_{i}^{2}} - (i-1)!*sumlimits_{j=1}^{i}frac{f_{j}}{(j-1)!} * g_{i-j} )

然后就可以分治FFT了；

很容易写错的地方是那个 ( 2^{C_{i}^{2}} )，因为是指数，所以应该对 ( mod-1 ) 取模，而不是 ( mod ) ！！

而且除以 ( 2 ) 也不能预处理 ( inv2 ) 了，因为 ( mod-1 ) 不是质数！！

代码如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int const xn=(1<<18),mod=1004535809;
int n,f[xn],g[xn],jc[xn],jcn[xn],a[xn],b[xn],rev[xn],inv2,in[xn];
int rd()
{
  int ret=0,f=1; char ch=getchar();
  while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();}
  while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
  return f?ret:-ret;
}
ll pw(ll a,int b,int p=mod)
{
  ll ret=1;
  for(;b;b>>=1,a=(a*a)%p)if(b&1)ret=(ret*a)%p;
  return ret;
}
int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;}
int C(int x){return ((ll)x*(x-1)/2)%(mod-1);}//mod-1!!!  //not inv2 -- !pri
void init()
{
  jc[0]=1;
  for(int i=1;i<=n;i++)jc[i]=(ll)jc[i-1]*i%mod;
  jcn[n]=pw(jc[n],mod-2);
  for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod;
  g[0]=0;
  for(int i=1;i<=n;i++)g[i]=(ll)pw(2,C(i))*jcn[i]%mod;
}
void ntt(int *a,int tp,int lim)
{
  for(int i=0;i<lim;i++)
    if(i<rev[i])swap(a[i],a[rev[i]]);
  for(int mid=1;mid<lim;mid<<=1)
    {
      int len=(mid<<1),wn=pw(3,tp==1?(mod-1)/len:(mod-1)-(mod-1)/len);
      for(int j=0;j<lim;j+=len)
    for(int k=0,w=1;k<mid;k++,w=(ll)w*wn%mod)
      {
        int x=a[j+k],y=(ll)w*a[j+mid+k]%mod;
        a[j+k]=upt(x+y); a[j+mid+k]=upt(x-y);
      }
    }
  if(tp==1)return; int inv=pw(lim,mod-2);
  for(int i=0;i<lim;i++)a[i]=(ll)a[i]*inv%mod;
}
void work(int l,int r)
{
  if(l==r){f[l]=upt((ll)pw(2,C(l))-(ll)jc[l-1]*f[l]%mod); return;}
  int len=r-l+1,mid=((l+r)>>1);
  work(l,mid);
  int lim=1,L=0;
  while(lim<=len)lim<<=1,L++;//
  for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(L-1)));
  for(int i=0;i<lim;i++)a[i]=b[i]=0;
  for(int i=l;i<=mid;i++)a[i-l]=(ll)f[i]*jcn[i-1]%mod;//jcn!!
  for(int i=0;i<=r-l;i++)b[i]=g[i];//
  ntt(a,1,lim); ntt(b,1,lim);
  for(int i=0;i<lim;i++)a[i]=(ll)a[i]*b[i]%mod;
  ntt(a,-1,lim);
  for(int i=mid+1;i<=r;i++)f[i]=upt(f[i]+a[i-l]);
  work(mid+1,r);
}
int main()
{
  n=rd(); init();
  work(1,n);
  printf("%d
",f[n]);
  return 0;
}

也可以用多项式求逆做：

设 ( g[i] = 2^{C_{i}^{2}} )，这里 ( g[0] = 1 )，( f[i] ) 定义同上；

得到 ( g[n] = sumlimits_{i=1}^{n} C_{n-1}^{i-1} * f[i] * g[n-i] )

拆 ( C_{n-1}^{i-1} )，得到

( frac{g[n]}{(n-1)!} = sumlimits_{i=1}^{n} frac{f[i]}{(i-1)!} * frac{g[n-i]}{(n-i)!} )

设 ( A[i] = frac{g[i]}{(i-1)!} ) ， ( B[i] = frac{f[i]}{(i-1)!} ) ， ( C[i] = frac{g[i]}{i!} )

注意这里 ( A[0] = 0 ) ， ( B[0] = 0 ) ， ( C[0] = 1 )

所以 ( A(x) = B(x) * C(x) )

即 ( B(x) = A(x) * C^{-1}(x) ( mod x^{n+1} ) )

多项式求逆即可；

别忘了外层的乘法还要处理一下 ( rev ) 数组。

代码如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int const xn=(1<<18),mod=1004535809;
int n,a[xn],b[xn],c[xn],t[xn],rev[xn],jc[xn],jcn[xn];
int rd()
{
  int ret=0,f=1; char ch=getchar();
  while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();}
  while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
  return f?ret:-ret;
}
ll pw(ll a,int b)
{
  ll ret=1;
  for(;b;b>>=1,a=(a*a)%mod)if(b&1)ret=(ret*a)%mod;
  return ret;
}
int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;}
int C(int x){return ((ll)x*(x-1)/2)%(mod-1);}
void init()
{
  jc[0]=1;
  for(int i=1;i<=n;i++)jc[i]=(ll)jc[i-1]*i%mod;
  jcn[n]=pw(jc[n],mod-2);
  for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod;
}
void ntt(int *a,int tp,int lim)
{
  for(int i=0;i<lim;i++)
    if(i<rev[i])swap(a[i],a[rev[i]]);
  for(int mid=1;mid<lim;mid<<=1)
    {
      int len=(mid<<1),wn=pw(3,tp==1?(mod-1)/len:(mod-1)-(mod-1)/len);
      for(int j=0;j<lim;j+=len)
    for(int k=0,w=1;k<mid;k++,w=(ll)w*wn%mod)
      {
        int x=a[j+k],y=(ll)w*a[j+mid+k]%mod;
        a[j+k]=upt(x+y); a[j+mid+k]=upt(x-y);
      }
    }
  if(tp==1)return; int inv=pw(lim,mod-2);
  for(int i=0;i<lim;i++)a[i]=(ll)a[i]*inv%mod;
}
void inv(int *a,int *b,int n)
{
  if(n==1){b[0]=pw(a[0],mod-2); return;}
  inv(a,b,(n+1)>>1);
  int lim=1,l=0;
  while(lim<n+n)lim<<=1,l++;
  for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));
  for(int i=0;i<n;i++)t[i]=a[i];
  for(int i=n;i<lim;i++)t[i]=0;
  ntt(t,1,lim); ntt(b,1,lim);
  for(int i=0;i<lim;i++)b[i]=upt((2-(ll)t[i]*b[i])%mod*b[i]%mod);
  ntt(b,-1,lim);
  for(int i=n;i<lim;i++)b[i]=0;
}
int main()
{
  n=rd(); init();
  int lim=1,l=0;
  while(lim<=n)lim<<=1,l++;
  for(int i=1;i<lim;i++)a[i]=(ll)pw(2,C(i))*jcn[i-1]%mod;
  for(int i=1;i<lim;i++)c[i]=(ll)pw(2,C(i))*jcn[i]%mod;
  c[0]=1;
  inv(c,b,n+1);
  for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));//!!!
  ntt(a,1,lim); ntt(b,1,lim);
  for(int i=0;i<lim;i++)b[i]=(ll)a[i]*b[i]%mod;
  ntt(b,-1,lim);
  printf("%lld
",(ll)b[n]*jc[n-1]%mod);
  return 0;
}

也可以用指数型生成函数，具体做法可以看这篇博客：https://blog.csdn.net/wzq_qwq/article/details/48435621

关于为什么 ( G(x) = e^{F(x)} )

因为从组合意义上来看，任意无向图实际上可以分为：由0个连通图组成(没有点)，由1个连通图组成，由2个连通图组成(这2个之间不连通)，由3个连通图组成...

所以 ( G(x) = 1 + frac{F(x)}{1!} + frac{F(x)^{2}}{2!} + frac{F(x)^{3}}{3!} + ... )

即 ( G(x) = e^{F(x)} )

而 ( G(x) ) 很好构造，求 ( F(x) = lnG(x) )

求 ( ln ) 的方法可以看这两篇博客：

https://blog.csdn.net/litble/article/details/81749788

https://blog.csdn.net/ezoixx118/article/details/81235586

直接用上公式...这题原来是多项式求 ( ln ) 的裸题；

别忘了最后乘上 ( n! )

代码如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int const xn=(1<<18),mod=1004535809;
int n,g[xn],dg[xn],ig[xn],t[xn],rev[xn],jc[xn],jcn[xn];
int rd()
{
  int ret=0,f=1; char ch=getchar();
  while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();}
  while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
  return f?ret:-ret;
}
ll pw(ll a,ll b)
{
  ll ret=1; b=b%(mod-1);
  for(;b;b>>=1,a=(a*a)%mod)if(b&1)ret=(ret*a)%mod;
  return ret;
}
int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;}
void init()
{
  jc[0]=1;
  for(int i=1;i<=n;i++)jc[i]=(ll)jc[i-1]*i%mod;
  jcn[n]=pw(jc[n],mod-2);
  for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod;
}
void ntt(int *a,int tp,int lim)
{
  for(int i=0;i<lim;i++)
    if(i<rev[i])swap(a[i],a[rev[i]]);
  for(int mid=1;mid<lim;mid<<=1)
    {
      int len=(mid<<1),wn=pw(3,tp==1?(mod-1)/len:(mod-1)-(mod-1)/len);
      for(int j=0;j<lim;j+=len)
    for(int k=0,w=1;k<mid;k++,w=(ll)w*wn%mod)
      {
        int x=a[j+k],y=(ll)w*a[j+mid+k]%mod;
        a[j+k]=upt(x+y); a[j+mid+k]=upt(x-y);
      }
    }
  if(tp==1)return; int inv=pw(lim,mod-2);
  for(int i=0;i<lim;i++)a[i]=(ll)a[i]*inv%mod;
}
void inv(int *a,int *b,int n)
{
  if(n==1){b[0]=pw(a[0],mod-2); return;}
  inv(a,b,(n+1)>>1);
  int lim=1,l=0;
  while(lim<n+n)lim<<=1,l++;
  for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));
  for(int i=0;i<n;i++)t[i]=a[i];
  for(int i=n;i<lim;i++)t[i]=0;
  ntt(t,1,lim); ntt(b,1,lim);
  for(int i=0;i<lim;i++)b[i]=upt((2-(ll)t[i]*b[i])%mod*b[i]%mod);
  ntt(b,-1,lim);
  for(int i=n;i<lim;i++)b[i]=0;
}
int main()
{
  n=rd(); init();
  for(int i=0;i<=n;i++)
    {
      if(i<2)g[i]=1;
      g[i]=(ll)pw(2,(ll)i*(i-1)/2)*jcn[i]%mod; 
    }
  for(int i=1;i<=n;i++)dg[i-1]=(ll)i*g[i]%mod;
  dg[n]=0;//
  inv(g,ig,n+1);
  int lim=1,l=0;
  while(lim<n+n)lim<<=1,l++;
  for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));
  ntt(ig,1,lim); ntt(dg,1,lim);
  for(int i=0;i<lim;i++)dg[i]=(ll)ig[i]*dg[i]%mod;
  ntt(dg,-1,lim);
  printf("%lld
",(ll)dg[n-1]*pw(n,mod-2)%mod*jc[n]%mod);
  return 0;
}