bzoj3456 城市规划(分治NTT)

本文介绍了一道经典算法题目“BZOJ3456城市规划”的解决思路,通过数学推导得到求解不同贸易路线方案数量的方法,并采用NTT算法进行高效计算。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

bzoj3456 城市规划

原题地址http://www.lydsy.com/JudgeOnline/problem.php?id=3456

题意:
阿狸的国家有n个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通. 为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案。

换句话说, 你需要求出n个点的简单(无重边无自环)无向连通图数目。

由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^ 21 + 1)即可。

数据范围
n <= 130000

题解:
gi g i i i 个点的图的数量,则gi=2n(n1)2,令 fi f i i i 个点的简单无向连通图数目,有:

fi=gij=1i1fj(i1j1)gij

对于减掉的那个式子:
j=1i1fj(i1j1)gij=j=1i1(i1)!(j1)!(ij)!fjgij=(i1)!j=1i1fj(j1)!gij(ij)! ∑ j = 1 i − 1 f j ∗ ( i − 1 j − 1 ) ∗ g i − j = ∑ j = 1 i − 1 ( i − 1 ) ! ( j − 1 ) ! ( i − j ) ! ∗ f j ∗ g i − j = ( i − 1 ) ! ∑ j = 1 i − 1 f j ( j − 1 ) ! ∗ g i − j ( i − j ) !
若令 Fi=fi(i1)! F i = f i ( i − 1 ) ! Gi=gij(ij)! G i = g i − j ( i − j ) !
fi=gi(i1)!j=1i1FjGij f i = g i − ( i − 1 ) ! ∑ j = 1 i − 1 F j ∗ G i − j ,如果让未算出的 Fn=F0=0 F n = F 0 = 0 ,且 G0=0 G 0 = 0 ,减掉的就完全转化为了标准卷积形式,可以使用NTT了:
fi=gi(i1)!j=0iFjGij f i = g i − ( i − 1 ) ! ∑ j = 0 i F j ∗ G i − j
由于计算出 fi f i 需要用到 f0...fi1 f 0 . . . f i − 1 ,所以用分治NTT即可。

说起来好像很简单的样子,然而调题调一天…

第一个错误:wi没有赋初值为1。
第二个错误:预处理g数组时因为指数是LL乘爆了。发现这个问题是因为 44824这组数据过不去。

要注意的地方大概是分治时计算[lf,mid]的f对[mid+1,rg]的贡献时,要想清楚哪一位具体对应哪一位。因为lf需要加上(rg-lf)才等于rg,所以用于相乘的g数组是[1,rg-lf]。

代码:

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#define LL long long
using namespace std;
const int mod=1004535809;//479*2^21+1 g=3
const int N=600050;
int G=3;
int n,fac[N],inv[N],f[N],a[N],g[N],b[N],R[N],iv[N],w[N],_w[N];
int modpow(int A,LL B)
{
    int ret=1; int base=A;
    for(;B;B>>=1)
    {
        if(B&1) ret=(1LL*ret*base)%mod;
        base=(1LL*base*base)%mod;
    }
    return ret;
}
void NTT(int *x,int flag,int len)
{
    for(int i=0;i<len;i++) if(i<R[i]) swap(x[i],x[R[i]]);
    for(int m=2;m<=len;m<<=1)
    {
        int l=m>>1; int wn,wi=1;    
        if(flag==1) wn=w[m]; else wn=_w[m];
        for(int j=0;j<len;j+=m)
        {
            wi=1;
            for(int i=0;i<l;i++)
            {
                int y=(1LL*wi*x[i+j+l])%mod;
                x[i+j+l]=(x[i+j]-y+mod)%mod;
                x[i+j]=(x[i+j]+y)%mod;
                wi=(1LL*wi*wn)%mod;
            }
        }

    }
    if(flag==-1) for(int i=0;i<len;i++) x[i]=(1LL*x[i]*iv[len])%mod;
}
void solve(int lf,int rg)
{
    if(lf==rg) {f[lf]=(g[lf]-(1LL*f[lf]*fac[lf-1])%mod+mod)%mod; return;}
    int mid=(lf+rg)>>1;
    solve(lf,mid);
    int len=1,p=0; for(;len<2*(rg-lf+1);len<<=1,p++);
    R[0]=0; for(int i=1;i<len;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(p-1));
    for(int i=0;i<=len;i++) a[i]=b[i]=0;
    for(int i=lf;i<=mid;i++) a[i-lf+1]=(1LL*f[i]*inv[i-1])%mod;
    for(int i=1;i<=rg-lf;i++) b[i]=(1LL*g[i]*inv[i])%mod;
    NTT(a,1,len); NTT(b,1,len);
    for(int i=0;i<len;i++) a[i]=(1LL*a[i]*b[i])%mod;
    NTT(a,-1,len);
    for(int i=mid+1;i<=rg;i++) f[i]=(f[i]+a[i-lf+1])%mod;
    solve(mid+1,rg);
}
int main()
{
    scanf("%d",&n); 
    inv[0]=fac[0]=inv[1]=fac[1]=1;
    for(int i=2;i<=4*n;i++) fac[i]=1LL*fac[i-1]*i%mod,inv[i]=modpow(fac[i],mod-2);
    g[0]=0; for(int i=1;i<=4*n;i++) g[i]=modpow(2,(1LL*i*(i-1)/2LL));
    for(int i=1,p=0;i<=4*n;i<<=1,p++)
    {w[i]=modpow(G,(mod-1)/i),_w[i]=modpow(w[i],mod-2),iv[i]=modpow(i,mod-2);}
    solve(1,n);
    printf("%d\n",f[n]);
    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值