jzoj3303. 【集训队互测2013】城市规划

Description

刚刚解决完电力网络的问题, 阿狸又被领导的任务给难住了.
刚才说过, 阿狸的国家有n 个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通.
为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.
好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出n 个点的简单(无重边无自环)无向连通图数目.
由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^21 + 1)即可.

Input

仅一行一个整数n(<=130000)

Output

仅一行一个整数, 为方案数mod 1004535809.

Sample Input

输入1:
3

输入2:
4

输入3:
100000

Sample Output

输出1:
4

输出2:
38

输出3:
829847355

Data Constraint

对于20%的数据, n <= 10
对于40%的数据, n <= 1000
对于60%的数据, n <= 30000
对于80%的数据, n <= 60000
对于100%的数据, n <= 130000

题解

题目疯狂暗示NTT。
然鹅比赛时不会,基本属于弃疗状态。
比赛完后疯狂学习FFT与NTT。
总算有个了解了。
不懂NTT的点这个:https://blog.csdn.net/HiChocolate/article/details/95023398

现在我们了解NTT是什么东东了。
回到这一题:
10%
直接暴力
30%
首先我们考虑DP。
直接DP求答案不是那么好搞,考虑容斥
f[i]f_{[i]}表示当前图大小为i的方案数。
方程长这样:
f[n]=2Cn2i=1n1f[i]Cn1i12Cni2f_{[n]}=2^{C_n^2}-sum_{i=1}^{n-1}f_{[i]}*C_{n-1}^{i-1}*2^{C_{n-i}^2}
为什么,请诸君一一听我细细道来:
前面的2Cn22^{C_n^2}表示当前无向图中所有的边,连接或不连接的方案数。
在sigma里面表示:我们钦定一个点在一个合法的联通图内,然后其他的点随意选择。
十分显然明了。
时间复杂度:O(n2)O(n^2)
话说某xzb强行把O(n2)O(n^2)的方法卡过了,%%%orz。
100%
我们看到上面的柿子。
考虑把那个丑陋的Cn1i1C_{n-1}^{i-1}给拆掉
得到:
f[n]=2Cn2i=1n1f[i](n1)!(i1)!(ni)!2Cni2f_{[n]}=2^{C_n^2}-sum_{i=1}^{n-1}f_{[i]}*frac{(n-1)!}{(i-1)!(n-i)!}*2^{C_{n-i}^2}
移项得:
f[n]=2Cn2(n1)!i=1n1f[i](i1)!2Cni2(ni)!f_{[n]}=2^{C_n^2}-(n-1)!*sum_{i=1}^{n-1}frac{f_{[i]}}{(i-1)!}*frac{2^{C_{n-i}^2}}{(n-i)!}
我们设A=f[i](i1)!frac{f_{[i]}}{(i-1)!}B=2Cni2(ni)!frac{2^{C_{n-i}^2}}{(n-i)!}
那么上面这个柿子可以表示成一个卷积的形式。
NTT打法好!
但是怎么套呢?
两个方法:
cdq分治
我们发现,当前的n都是有前面的i得到的。
那么考虑分治。
每次分治一个区间[l,mid]把这个区间的代价加到[mid+1,r]里面去。
怎么加?
将[l,mid]的A提出来,与一个长度为r-l的B跑一边NTT,得到一个新的C,加入[mid+1,r]里面。
对于一个单点,直接求解答案即可。
思路清晰,实现恶心。
多项式求逆
还是画柿子。
稍微变一下:
2Cn2(n1)!=i=1nf[i](i1)!2Cni2(ni)!frac{2^{C_n^2}}{(n-1)!}=sum_{i=1}^{n}frac{f_{[i]}}{(i-1)!}*frac{2^{C_{n-i}^2}}{(n-i)!}
这个是个更加明显的卷积形式。
A=2Cn2(n1)!B=f[i](i1)!C=2Cni2(ni)!A=frac{2^{C_n^2}}{(n-1)!}B=frac{f_{[i]}}{(i-1)!}C=frac{2^{C_{n-i}^2}}{(n-i)!}
只不过里面要求的B在卷积里面。
怎么办?
我们换一下即可知道,我们求出C的逆序列即可再次利用NTT。
这叫做多项式求逆。
这里有个blog,自己学。(我不想打)
https://www.cnblogs.com/yoyoball/p/8724115.html

代码

cdq分治的(跑得贼慢)

{$inline on}
uses math;
type
        arr=array[-10..262144] of int64;
var
        i,j,k,l,n,m,lim:longint;
        rev:int64;
        mo:int64=1004535809;
        f,g,w,a,b,jc,ny:arr;
function qsm(a,b:int64):int64;
var
        t,y:int64;
begin
        t:=1;
        y:=a;
        while b>0 do
        begin
                if b and 1=1 then t:=t*y mod mo;
                y:=y*y mod mo;
                b:=b div 2;
        end;
        exit(t);
end;
procedure dft(var a:arr;n:longint);inline;
var
        i,j,k,now:longint;
        x,y:int64;
        r:arr;
begin
        k:=n-1;j:=0;
        while k>0 do
        begin
                k:=k div 2;
                inc(j);
        end;
        r[0]:=0;
        for i:=0 to n-1 do
        begin
                r[i]:=(r[i div 2] div 2) or ((i and 1) shl (j-1));
        end;
        for i:=0 to n-1 do
        begin
                if i<r[i] then
                begin
                        j:=a[i];a[i]:=a[r[i]];a[r[i]]:=j;
                end;
        end;
        now:=2;
        while now<=n do
        begin
                i:=now div 2;
                for j:=0 to i-1 do
                begin
                        k:=j;
                        while k<n do
                        begin
                                x:=a[k];y:=a[k+i]*w[j*(n div now)] mod mo;
                                a[k]:=(x+y) mod mo;
                                a[k+i]:=(x-y+mo) mod mo;
                                k:=k+now;
                        end;
                end;
                now:=now*2;
        end;
end;
procedure fft(var a,b:arr;n:longint);
var
        i:longint;
        j,k:int64;
begin
        rev:=qsm(3,(mo-1) div n);
        w[0]:=1;
        for i:=1 to n do w[i]:=w[i-1]*rev mod mo;
        dft(a,n);
        dft(b,n);
        for i:=0 to n-1 do a[i]:=a[i]*b[i] mod mo;
        for i:=0 to n div 2 do
        begin
                j:=w[i];w[i]:=w[n-i];w[n-i]:=j;
        end;
        dft(a,n);
        k:=qsm(n,mo-2);
        for i:=0 to n-1 do
        begin
                a[i]:=a[i]*k mod mo;
        end;
end;
procedure solve(l,r:longint);
var
        i,j,k,mid:longint;

begin
        if l=r then
        begin
                if r=44824 then
                mid:=mid;
                f[l]:=(g[l]-f[l]*jc[l-1] mod mo+mo) mod mo;
        end
        else
        begin
                mid:=(l+r) div 2;
                solve(l,mid);
                if mid-l+1>r-l then j:=mid-l+1
                else j:=r-l;
                k:=1;
                while k<=j do k:=k*2;
                for i:=0 to 2*k do
                begin
                        a[i]:=0;
                        b[i]:=0;
                end;
                for i:=l to mid do a[i-l+1]:=f[i]*ny[i-1] mod mo;
                for i:=1 to r-l do b[i]:=g[i]*ny[i] mod mo;
                j:=max((mid-l+1)*2,r-l+1);
                fft(a,b,j*2);
                for i:=mid+1 to r do
                begin
                        f[i]:=(f[i]+a[i-l+1]) mod mo;
                end;
                solve(mid+1,r);
        end;
end;
begin
        //assign(input,'4data.in');reset(input);
        g[0]:=1;
        for i:=1 to 262144 do g[i]:=qsm(2,((int64(i-1)*i div 2) mod (mo-1)) mod mo);
        jc[0]:=1;
        for i:=1 to 262144 do jc[i]:=jc[i-1]*i mod mo;
        ny[262144]:=qsm(jc[262144],mo-2);
        for i:=262143 downto 0 do ny[i]:=ny[i+1]*(i+1) mod mo;
        readln(n);
        lim:=1;
        while lim<n+m do lim:=lim*2;
        solve(1,lim);
        writeln(f[n]);
end.
end.



原文地址:https://www.cnblogs.com/RainbowCrown/p/11148345.html