HDU 5896 Running King(dp+NTT+CDQ分治)

本文介绍了一种计算特定数量顶点的无向有环图数量的方法,并提供了一个使用快速傅里叶变换(FFT)优化的算法实现。该算法通过递归分解问题规模并利用卷积技巧来减少计算复杂度。

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

Description

nn个点的无向有环图个数

Input

第一行一整数T表示用例组数,每组用例输入一整数n(1T10,1n2105)n(1≤T≤10,1≤n≤2⋅105)

Output

输出nn个点的无向有环图个数,结果模1004535809

Sample Input

3
2
3
4

Sample Output

0
1
26

Solution

nn个点的无向有环图个数=2n(n1)2n个点的无向森林个数

f(n)f(n)nn个点的无向森林个数,g(n)nn个点的无向树的个数,则g(n)=nn2(每棵nn个节点无向树唯一对应一个长度为n2PurferPurfer序列)

枚举第一个点所在树的节点数有转移f(n)=g(n)+i=1n1Ci1n1g(i)f(ni)f(n)=g(n)+∑i=1n−1Cn−1i−1g(i)f(n−i),令f(0)=1f(0)=1,则有f(n)(n1)!=i=1nii1i!f(ni)(ni)!f(n)(n−1)!=∑i=1nii−1i!f(n−i)(n−i)!CDQCDQ分治即可

假设当前求解f(l),...,f(r)f(l),...,f(r),而f(l),...,f(mid)f(l),...,f(mid)已经求出,考虑f(l),...,f(mid)f(l),...,f(mid)f(mid+1),...,f(r)f(mid+1),...,f(r)的贡献,令x(i)=(i+1)i(i+1)!,y(i)=f(i+l)(i+l)!,i=0,...,midlx(i)=(i+1)i(i+1)!,y(i)=f(i+l)(i+l)!,i=0,...,mid−l,对mid+1trmid+1≤t≤rf(l),...,f(mid)f(l),...,f(mid)f(t)(t1)!f(t)(t−1)!的贡献为h(t)=i=lmidx(ti1)f(i)i!=i=0tl1x(tl1i)y(i)=(xy)(tl1)h(t)=∑i=lmidx(t−i−1)f(i)i!=∑i=0t−l−1x(t−l−1−i)y(i)=(x∗y)(t−l−1),其中xyx∗y表示x,yx,y两序列的卷积,NTTNTT即可快速求卷积

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long ll;
#define maxn 222222
typedef long long ll;
const ll mod=1004535809ll;
const ll g=3;
ll mod_pow(ll a,ll b,ll p)
{
    a%=p;
    ll ans=1;
    while(b)
    {
        if(b&1)ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
void change(ll *x,int len) 
{
    for(int i=1,j=len/2;i<len-1;i++)
    {
        if(i<j)swap(x[i],x[j]);
        int k=len/2;
        while(j>=k)
        {
            j-=k;
            k/=2;
        }
        if(j<k)j+=k;
    }
}
void ntt(ll *x,int len,int sta) 
{
    change(x,len);
    for(int m=2;m<=len;m<<=1) 
    {
        ll Wn=mod_pow(g,(mod-1)/m,mod);
        if(sta==-1)Wn=mod_pow(Wn,mod-2,mod);
        for(int i=0;i<len;i+=m) 
        {
            ll W=1;
            for(int j=i;j<i+m/2;j++) 
            {
                ll x1=x[j],x2=W*x[j+m/2]%mod;
                x[j]=(x1+x2)%mod,x[j+m/2]=(x1-x2+mod)%mod;
                W=W*Wn%mod;
            }
        }
    }
    if(sta==-1)
    {
        int temp=mod_pow(len,mod-2,mod);
        for(int i=0;i<len;i++)
            x[i]=x[i]*temp%mod;
    }
}
ll mod_pow(ll a,ll b)
{
    ll ans=1ll;
    while(b)
    {
        if(b&1)ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
} 
ll n,a[maxn],dp[maxn];
ll x[2*maxn],y[2*maxn];
ll f[maxn],h[maxn],inv[maxn];
void init()
{
    f[0]=1,inv[0]=1;
    for(int i=1;i<maxn;i++)f[i]=f[i-1]*i%mod,inv[i]=mod_pow(f[i],mod-2);
    for(int i=0;i<maxn;i++)h[i]=mod_pow(i+1,i);
}
void deal(int l,int r)
{
    if(l==r)return ;
    int mid=(l+r)>>1;
    deal(l,mid);
    int len=1;
    while(len<=r-l+1)len<<=1;
    for(int i=0;i<len;i++)
    {
        x[i]=a[i];
        if(i+l<=mid)y[i]=dp[i+l]*inv[i+l]%mod;
        else y[i]=0;
    }
    ntt(x,len,1),ntt(y,len,1);
    for(int i=0;i<len;i++)x[i]=x[i]*y[i]%mod;
    ntt(x,len,-1);
    for(int i=mid+1;i<=r;i++)
        dp[i]+=x[i-l-1]*f[i-1]%mod,dp[i]%=mod;
    deal(mid+1,r);
}
int main()
{
    int T;
    scanf("%d",&T);
    init();
    while(T--)
    {
        scanf("%d",&n);
        for(int i=0;i<=n;i++)a[i]=h[i]*inv[i+1]%mod;
        memset(dp,0,sizeof(dp));
        dp[0]=1;
        deal(0,n);
        ll ans=mod_pow(2,1ll*n*(n-1)/2)-dp[n];
        printf("%I64d\n",(ans%mod+mod)%mod);
    }
    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值