bzoj 3456 城市规划 - 图计数 - NTT - 多项式求逆 - NTT学习笔记 - 多项式求逆学习笔记

64 篇文章 0 订阅
35 篇文章 0 订阅

设f[x]表示x个点的连通图,h[x]=2^{C(n,2)}为任意图的数量。
H n = ∑ i = 1 n ( n − 1 i − 1 ) F i ×   H n − i H_n=\sum_{i=1}^n \binom{n-1}{i-1}F_i\times\ H_{n-i} Hn=i=1n(i1n1)Fi× Hni
H n ( n − 1 ) ! = ∑ i = 1 n F i ( i − 1 ) ! ×   H n − i ( n − i ) ! \frac{H_n}{(n-1)!}=\sum_{i=1}^n \frac{F_i}{(i-1)!}\times\ \frac{H_{n-i}}{(n-i)!} (n1)!Hn=i=1n(i1)!Fi× (ni)!Hni
A ( x ) = ∑ i = 0 n H n ( n − 1 ) ! x i , F ( x ) = ∑ i = 0 n F i ( i − 1 ) ! x i , B ( x ) = ∑ i = 0 n H n n ! x i A(x)=\sum_{i=0}^n\frac{H_n}{(n-1)!}x^i, F(x)=\sum_{i=0}^n\frac{F_i}{(i-1)!}x^i,B(x)=\sum_{i=0}^n\frac{H_n}{n!}x^i A(x)=i=0n(n1)!Hnxi,F(x)=i=0n(i1)!Fixi,B(x)=i=0nn!Hnxi
那么 A = F × B ,   F = A × B − 1 A=F\times B,\ F=A\times B^{-1} A=F×B, F=A×B1。多项式求逆即可。
代码:

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define N 1200010
#define lint long long
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
int A[N],B[N],C[N],D[N],a[N],b[N],c[N],x[N],r[N],fac[N],facinv[N];
inline int show(int *A,int n)
{
    for(int i=0;i<n;i++) cerr<<A[i]<<" ";cerr ln;
    return 0;
}
inline int fast_pow(int x,lint k,int p)
{
    int ans=1;
    while(k) ((k&1ll)?ans=(lint)ans*x%p:0),x=(lint)x*x%p,k>>=1ll;
    return ans;
}
inline int prelude(int n,int p)
{
    for(int i=fac[0]=1;i<=n;i++) fac[i]=(lint)fac[i-1]*i%p;
    facinv[n]=fast_pow(fac[n],p-2,p);
    for(int i=n-1;i>=0;i--) facinv[i]=(lint)facinv[i+1]*(i+1)%p;
    return 0;
}
inline int NTT(int *a,int n,int sgn,int p=1004535809,int g=3)
{
    int L=0;for(int i=1;i<n;i<<=1,L++);
    for(int i=1;i<=n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
    for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=2;i<=n;i<<=1)
    {
        int wn=fast_pow(g,(sgn>0)?((p-1)/i):(p-1-(p-1)/i),p);
        for(int j=0,t=i>>1;j<n;j+=i)
            for(int k=0,w=1;k<t;k++,w=(lint)w*wn%p)
            {
                int x=a[j+k],y=(lint)w*a[j+k+t]%p;
                a[j+k]=(x+y)%p,a[j+k+t]=(x-y+p)%p;
            }
    }
    return 0;
}
inline int solve(int *A,int *B,int *c,int n,int p)
{
    for(int i=0;i<n;i++) a[i]=A[i];
    for(int i=0;i<n;i++) b[i]=B[i];
    NTT(a,n,1,p,3),NTT(b,n,1,p,3);
    for(int i=0;i<n;i++) c[i]=(lint)a[i]*b[i]%p;
    NTT(c,n,-1,p,3);int ninv=fast_pow(n,p-2,p);
    for(int i=0;i<n;i++) c[i]=(lint)c[i]*ninv%p;
    return 0;
}
inline int Polynomial_Inv(int *A,int *B,int n,int p)
{
    B[0]=1;
    for(int i=2;i<=n;i<<=1)
    {
        for(int j=0;j<i;j++) x[j]=2ll*B[j]%p;solve(B,B,B,i,p);
        for(int j=0;j<i;j++) c[j]=A[j];solve(c,B,B,i<<1,p);
        for(int j=0;j<i;j++) B[j]=(x[j]-B[j]+p)%p;
        for(int j=i;j<(i<<1);j++) B[j]=0;
    }
    return 0;
}
inline int Init(int n,int p)
{
    prelude(n,p);
    for(int i=0;i<n;i++) B[i]=(lint)fast_pow(2,(lint)i*(i-1)/2,p)*facinv[i]%p;
    for(int i=1;i<n;i++) C[i]=(lint)fast_pow(2,(lint)i*(i-1)/2,p)*facinv[i-1]%p;
    return 0;
}
inline int calc(int n,int p=1004535809)
{
    int ans_pos=n;for(n=1;n<=ans_pos;n<<=1);
    Init(n,p),Polynomial_Inv(B,D,n,p),solve(D,C,A,n<<1,p);
    return (lint)A[ans_pos]*fac[ans_pos-1]%p;
}
int main()
{
    int n;scanf("%d",&n);
    return !printf("%d\n",calc(n));
}

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值