隐藏
「bzoj3456」城市规划 - 动态规划+生成函数+多项式求逆+NTT | Bill Yang's Blog

路终会有尽头,但视野总能看到更远的地方。

0%

「bzoj3456」城市规划 - 动态规划+生成函数+多项式求逆+NTT

题目大意

    求出有$n(n\le130000)$个点的有标号简单连通无向图的个数,答案$\bmod\,1004535809(479\times 221+1)$,是个质数。


题目分析

Miskcoo的题解讲得很详细了,不需要再赘述一遍。
实际上Miskoo手推了一遍指数型生成函数的式子。
注意我的多项式求逆与Miskoo的略有不同,我保证了$n$是$2$的整数幂,因此没有向上取整操作。
多项式求逆学习笔记

因为素数特殊的原因,故使用NTT而不使用FFT。

本题还有使用分治FFT的做法,不过较难调试,留一个坑。


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

typedef long long LL;

inline const int Get_Int() {
int num=0,bj=1;
char x=getchar();
while(x<'0'||x>'9') {
if(x=='-')bj=-1;
x=getchar();
}
while(x>='0'&&x<='9') {
num=num*10+x-'0';
x=getchar();
}
return num*bj;
}

const int maxn=262144+5;
const LL mod=1004535809,g=3;

LL Quick_Pow(LL a,LL b) {
LL sum=1;
for(; b; b>>=1,a=a*a%mod)if(b&1)sum=sum*a%mod;
return sum;
}

LL inv(LL x) {
return Quick_Pow(x,mod-2);
}

struct NumberTheoreticTransform {
int n;
LL omega[maxn],iomega[maxn];

void init(int n) {
this->n=n;
int x=Quick_Pow(g,(mod-1)/n);
omega[0]=iomega[0]=1;
for(int i=1; i<n; i++) {
omega[i]=omega[i-1]*x%mod;
iomega[i]=inv(omega[i]);
}
}

void transform(LL* a,LL* omega) {
int k=log2(n);
for(int i=0; i<n; i++) {
int t=0;
for(int j=0; j<k; j++)if(i&(1<<j))t|=(1<<(k-j-1));
if(i<t)swap(a[i],a[t]);
}
for(int len=2; len<=n; len*=2) {
int mid=len>>1;
for(LL* p=a; p!=a+n; p+=len)
for(int i=0; i<mid; i++) {
LL t=omega[n/len*i]*p[mid+i]%mod;
p[mid+i]=(p[i]-t+mod)%mod;
p[i]=(p[i]+t)%mod;
}
}
}

void dft(LL* a) {
transform(a,omega);
}

void idft(LL* a) {
transform(a,iomega);
LL x=inv(n);
for(int i=0; i<n; i++)a[i]=a[i]*x%mod;
}
} ntt;

void polynomial_inverse(const LL* a,const int n,LL* b) {
if(n==1) {
b[0]=inv(a[0]);
return;
}
polynomial_inverse(a,n>>1,b);
int p=n<<1;
static LL x[maxn];
copy(a,a+n,x),fill(x+n,x+p,0);
ntt.init(p),ntt.dft(x),ntt.dft(b);
for(int i=0; i<p; i++)b[i]=b[i]*((2-x[i]*b[i]%mod+mod)%mod)%mod;
ntt.idft(b),fill(b+n,b+p,0);
}

LL n,fac[maxn],invi[maxn],invf[maxn],g1[maxn],G[maxn],C[maxn],_G[maxn];

int main() {
n=Get_Int();
fac[0]=invf[0]=invi[1]=1;
for(int i=1; i<=n; i++) {
fac[i]=fac[i-1]*i%mod;
if(i!=1)invi[i]=(mod-mod/i)*invi[mod%i]%mod;
invf[i]=invf[i-1]*invi[i]%mod;
}
g1[0]=g1[1]=1;
for(int i=2; i<=n; i++)g1[i]=Quick_Pow(2,(LL)i*(i-1)/2%(mod-1));
for(int i=0; i<=n; i++)G[i]=g1[i]*invf[i]%mod;
for(int i=1; i<=n; i++)C[i]=g1[i]*invf[i-1]%mod;
int p=1;
while(p<n+1)p<<=1;
polynomial_inverse(G,p,_G);
fill(_G+n+1,_G+p,0);
p<<=1;
ntt.init(p);
ntt.dft(_G);
ntt.dft(C);
for(int i=0; i<=p; i++)_G[i]=_G[i]*C[i]%mod;
ntt.idft(_G);
printf("%lld\n",_G[n]*fac[n-1]%mod);
return 0;
}
姥爷们赏瓶冰阔落吧~