bzoj 3456 城市规划——分治FFT / 多项式求逆 / 多项式求ln

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

分治FFT：

　　设 dp[ i ] 表示 i 个点时连通的方案数。

　　考虑算补集：连通的方案数 == 随便连方案数 - 不连通方案数

　　不连通方案数就和很久之前做过的“地震后的幻想乡”一样，枚举一个连通的点集，其中需要一直包含一个“划分点”保证不重复；其余部分随便连。注意还有从 i 个点里选 j 个点作为连通点集的那个组合数。

　　( dp[i]=2^{C^{2}_{i}} - sumlimits^{i-1}_{j=1} dp[j]*C^{j-1}_{i-1}*2^{C^{2}_{i-j}} )

　　( dp[i]=2^{C^{2}_{i}} - (i-1)!sumlimits^{i-1}_{j=1} ( dp[j]*frac{1}{(j-1)!} )( 2^{C^{2}_{i-j}}*frac{1}{(i-j)!} ) )

　　就可以分治FFT啦！

　　一开始还记得指数是对 mod-1 取模，后来就忘了，调了好久……注意 mod-1 不是质数，且是偶数，所以2没有逆元，算 C 的时候直接 /2 就行了。

　　在 L==R 的地方把 f [ ] 弄好。虽然也可以给 f 赋上 ( C^{2}_{i} ) 的初值，然后卷积的时候每次 -=( (i-1)!*a[j] )，就不用在 L==R 的时候特殊判断了，但那样可能因为总要乘 ( (i-1)! )，会变慢。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int N=130005,M=N<<1,mod=1004535809,m2=mod-1;//N<<1
int n,len,r[M],f[M],g[M],a[M],b[M],jc[N],jcn[N],inv2;
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;
}
void upd(int &x){x>=mod?x-=mod:0;}
int pw(int x,int k,int md=mod)
{int ret=1;while(k){if(k&1)ret=(ll)ret*x%md;x=(ll)x*x%md;k>>=1;}return ret;}
int C(int n)///////%mod-1////
{return (ll)n*(n-1)/2%m2;}// /2
void init()
{
  //////  inv2=pw(2,m2-2,m2);//////m2!!!! m2 is not prime!!!!!
  jc[1]=1;for(int i=2;i<n;i++)jc[i]=(ll)jc[i-1]*i%mod;
  jcn[n-1]=pw(jc[n-1],mod-2);
  for(int i=n-2;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod;//>=0
  for(int i=1;i<n;i++)g[i]=(ll)pw(2,C(i))*jcn[i]%mod;
}
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=pw( 3,fx?(mod-1)-(mod-1)/R:(mod-1)/R );
      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]=x+y;  upd(a[i+j]);
        a[i+m+j]=x+mod-y;  upd(a[i+m+j]);
      }
    }
  if(!fx)return; int inv=pw(len,mod-2);
  for(int i=0;i<len;i++)a[i]=(ll)a[i]*inv%mod;
}
void solve(int L,int R)
{
  if(L==R){f[L]=(pw(2,C(L))-(ll)jc[L-1]*f[L])%mod+mod;upd(f[L]);return;}
  int mid=L+R>>1;  solve(L,mid);
  int d=R-L,i,j;
  for(len=1;len<d;len<<=1);
  for(i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)?len>>1:0);

  for(i=0,j=L;j<=mid;i++,j++)a[i]=(ll)f[j]*jcn[j-1]%mod;//jcn
  for(;i<len;i++)a[i]=0;
  for(i=0,j=R-L;i<=j;i++)b[i]=g[i+1];
  for(;i<len;i++)b[i]=0;
  ntt(a,0); ntt(b,0);
  for(i=0;i<len;i++)a[i]=(ll)a[i]*b[i]%mod;
  ntt(a,1);
  for(int i=mid+1,j=i-L-1;i<=R;i++,j++)f[i]+=a[j],upd(f[i]);
  solve(mid+1,R);
}
int main()
{
  n=rdn(); init();
  solve(1,n);
  printf("%d
",f[n]);
  return 0;
}

 多项式求逆：

　　http://blog.miskcoo.com/2015/05/bzoj-3456

　　应该只有G(x)有0次项。在阶乘的地方看，0!应该等于1才对（预处理组合数时就是这样，想一想，( C^{n}_{n} )就是( frac{n!}{n!*0!} )）。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int N=1e5+5,M=N<<2,mod=1004535809;
int n,a[M],c[M],cn[M],A[M],jcn[M],len,r[M];
void upd(int &x){x>=mod?x-=mod:0;}
int C(int n){return (ll)n*(n-1)/2%(mod-1);}
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;}
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=pw( 3,fx?(mod-1)-(mod-1)/R:(mod-1)/R );
      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]=x+y; upd(a[i+j]);
        a[i+m+j]=x+mod-y; upd(a[i+m+j]);
      }
    }
  if(!fx)return ; int inv=pw( len,mod-2 );
  for(int i=0;i<len;i++)a[i]=(ll)a[i]*inv%mod;
}
void getinv(int n,int *a,int *b)
{
  if(n==1){b[0]=1;return;}
  getinv(n+1>>1,a,b);
  int i;
  for(len=1;len<n<<1;len<<=1);
  for(i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)?len>>1:0);
  for(i=0;i<n;i++)A[i]=a[i];  for(;i<len;i++)A[i]=0;
  ntt(A,0);  ntt(b,0);
  for(i=0;i<len;i++)b[i]=(2-(ll)A[i]*b[i]%mod)*b[i]%mod+mod,upd(b[i]);
  ntt(b,1);
  for(int i=n;i<len;i++)b[i]=0;
}
int main()
{
  scanf("%d",&n);
  for(int i=1;i<=n;i++)a[i]=pw(2,C(i));
  jcn[0]=1;for(int i=1;i<=n;i++)jcn[i]=(ll)jcn[i-1]*i%mod;
  jcn[n]=pw(jcn[n],mod-2);for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod;
  c[0]=1; for(int i=1;i<=n;i++)c[i]=(ll)a[i]*jcn[i]%mod;
  for(int i=1;i<=n;i++)a[i]=(ll)a[i]*jcn[i-1]%mod;
  getinv(n+1,c,cn);
  for(len=1;len<=n<<1;len<<=1);
  for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)?len>>1:0);
  ntt(a,0); ntt(cn,0);
  for(int i=0;i<len;i++)a[i]=(ll)a[i]*cn[i]%mod;
  ntt(a,1);
  printf("%lld
",(ll)a[n]*pw(jcn[n-1],mod-2)%mod);
  return 0;
}

多项式求 ln ：

关于指数型生成函数：https://max.book118.com/html/2015/1111/29175835.shtm

想求两种物品一共取 i 个的排列方案，设取第一种物品的方案是 f[ ] ，第二种物品的方案是 g[ ] ，那么就是想求

　　( sumlimits_{j=0}^{i}C_{i}^{j}f[i]*g[i-j] )

令 F(x) 是 f[ ] 的指数型生成函数、G(x) 是 g[ ] 的生成函数的话，要求的就恰好是 F(x) * G(x) 的第 i 次项系数乘上 i! 。

关于本题：https://www.cnblogs.com/ivorysi/p/9756961.html

至于为什么式子 ( G(x)=sumlimits_{i=1}^{infty}frac{F(i)}{i!} ) 里要除以 i 的阶乘，可以从 n=2 并希望有两个连通块的情况考虑：

令 ( g[ i ] = 2^{C_{n}^{2}} ) 表示 i 个点的无向图个数，这个算出来是有标号的；令 f[ i ] 表示 i 个点有标号无向连通图个数。

( F^{2}(x) = sumlimits_{i=0}^{n}frac{ sumlimits_{j=0}^{i}C_{i}^{j}f[j]*f[i-j] }{i!}x^i )

直接这样的话，就有 ( g[2]=C_{2}^{0}*f[0]*f[2]+C_{2}^{1}*f[1]*f[1]+C_{2}^{2}f[2]*f[0] )

可以发现，因为有了那个 C 的系数，所以即使是卷积的时候只会算一遍的 j == i-j 的情况，也会有重复；即那个博客里说的本质一样的情况。

关于求ln：http://www.cnblogs.com/zhoushuyu/p/8763215.html

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int N=(1<<19)+5,mod=1004535809;
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,len,r[N],f[N],g[N],A[N],B[N];
int jc[N],jcn[N],inv[N];
void init(int n)
{
  jc[0]=jcn[0]=inv[0]=1;
  jc[1]=jcn[1]=inv[1]=1;
  for(int i=2;i<=n;i++)
    {
      jc[i]=(ll)jc[i-1]*i%mod;
      inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod;
      jcn[i]=(ll)jcn[i-1]*inv[i]%mod;
    }
}
void ntt(int *a,bool fx,bool deb=0)
{
  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=pw( 3,fx?(mod-1)-(mod-1)/R:(mod-1)/R );
      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 inv=pw(len,mod-2);
  for(int i=0;i<len;i++)a[i]=(ll)a[i]*inv%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]=b[n-1]=0;
}
void get_jf(int n,int *a)
{
  for(int i=n-1;i;i--)a[i]=(ll)a[i-1]*inv[i]%mod;//i--
  a[0]=0;
}
void get_inv(int n,int *a,int *b)
{
  if(n==1){b[0]=pw(a[0],mod-2);return;}
  get_inv(n>>1,a,b);
  len=n<<1;
  for(int i=0,j=len>>1;i<len;i++)
    r[i]=(r[i>>1]>>1)+((i&1)?j:0);
  for(int i=0;i<n;i++)A[i]=a[i],A[i+n]=0;
  ntt(A,0); ntt(b,0);
  for(int i=0;i<len;i++)
    b[i]=(ll)b[i]*(2-(ll)A[i]*b[i]%mod+mod)%mod;
  ntt(b,1); for(int i=n;i<len;i++)b[i]=0;
}
void get_ln(int n,int *a,int *b)
{
  get_inv(n,a,B); get_dao(n,a,A);//get_inv will use A
  len=n<<1; for(int i=n;i<len;i++)A[i]=0;//
  for(int i=0,j=len>>1;i<len;i++)
    r[i]=(r[i>>1]>>1)+((i&1)?j:0);
  ntt(A,0,1); ntt(B,0,1);
  for(int i=0;i<len;i++)b[i]=(ll)A[i]*B[i]%mod;
  ntt(b,1);
  get_jf(n,b);
  for(int i=n;i<len;i++)b[i]=0;//
}
int main()
{
  scanf("%d",&n);int l;for(l=1;l<=n;l<<=1);//l>n
  init(l);
  g[0]=1;///
  for(int i=1;i<=n;i++)
    g[i]=(ll)pw(2,(ll)i*(i-1)/2%(mod-1))*jcn[i]%mod;
  get_ln(l,g,f); printf("%lld
",(ll)f[n]*jc[n]%mod);
  return 0;
}