[BZOJ 3456]城市规划
aji多项式求逆毁我青春,,,
设fi表示i个点的有标号无向联通图,考虑所有可能的图(记Fi为i个点的有标号无向图,显然F_i = 2^{\binom{i}{2}})和f的关系(使用图计数的经典套路:枚举1所在的联通块大小):
F_n = \sum_{i = 1}^n \binom{n-1}{i-1} f_i F_{n - i}
看起来事卷积?但是这个卷积没有办法直接用FFT/NTT求(当然分离一下项啊,移下项就可以分治NTT力)。
考虑进一步化简柿子。完全展开后会发现右边有个非常碍眼的(n-1)!,所以两边除一下:
\frac{F_n}{(n-1)!} =\sum_{i = 1}^n \frac{f_i}{(i-1)!}\cdot \frac{F_{n - i}}{(n - i)!}
然后这个问题就很毒瘤了:我们要求的答案多项式(除上那个阶乘)和一个已知多项式做卷积,可以得到另一个已知多项式……这样就需要多项式除法了,于是乎多项式逆元派上了用场。
代码:
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 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 | #include <cstdio> #include <cstring> #include <cstdlib> #include <cctype> #include <algorithm> #include <utility> typedef long long ll; const int maxn = 1020010; const ll ha = 1004535809LL; const ll bs = 3LL; ll pow_mod(ll a, ll b) { ll ans = 1LL, res = a % ha; while (b) { if (1LL & b) ans = (ans * res) % ha; res = (res * res) % ha; b >>= 1; } return ans; } ll inv(ll x) { return pow_mod(x, ha - 2LL); } int flip( int bi, int x) { int ans = 0; for ( int i = 0; i < bi; i ++) { if ((1 << i) & x) { ans += (1 << (bi - i - 1)); } } return ans; } void NTT(ll *A, int bi, bool flag = false ) { int n = 1 << bi; for ( int i = 0; i < n; i ++) { int v = flip(bi, i); if (v < i) std::swap(A[v], A[i]); } for ( int L = 1; L < n; L <<= 1) { ll xi = pow_mod(3LL, (ha - 1LL) / (ll(L << 1))); if (flag) xi = inv(xi); for ( int i = 0; i < n; i += (L << 1)) { ll w = 1LL; for ( int j = i; j < i + L; j ++) { ll v1 = A[j], v2 = A[j + L]; A[j] = (v1 + (w * v2) % ha) % ha; A[j + L] = (v1 - (w * v2) % ha + ha) % ha; w = (w * xi) % ha; } } } } void poly_mul(ll *A, ll *B, int bi, ll *C) { static ll T1[maxn], T2[maxn]; int n = (1 << bi); std::copy(A, A + n, T1); std::copy(B, B + n, T2); NTT(T1, bi); NTT(T2, bi); #ifdef LOCAL puts ( "poly_mul :" ); for ( int i = 0; i < (n); i ++) { printf ( "%lld " , A[i]); } puts ( "" ); for ( int i = 0; i < (n); i ++) { printf ( "%lld " , B[i]); } puts ( "" ); for ( int i = 0; i < (n); i ++) { printf ( "%lld " , T1[i]); } puts ( "" ); for ( int i = 0; i < (n); i ++) { printf ( "%lld " , T2[i]); } puts ( "" ); #endif for ( int i = 0; i < n; i ++) { T1[i] = (T1[i] * T2[i]) % ha; } #ifdef LOCAL for ( int i = 0; i < (n); i ++) { printf ( "%lld " , T1[i]); } puts ( "" ); #endif NTT(T1, bi, true ); ll inv_n = inv(n); for ( int i = 0; i < n; i ++) { T1[i] = (T1[i] * inv_n) % ha; } std::copy(T1, T1 + n, C); #ifdef LOCAL for ( int i = 0; i < (n); i ++) { printf ( "%lld " , C[i]); } puts ( "" ); #endif } void poly_inv( int mod, ll *B, ll *BB) { if (mod == 1) { BB[0] = inv(B[0]); } else { poly_inv((mod + 1) >> 1, B, BB); int bi = 0, sz = 1; while (sz <= ((mod * 2) + 1)) { bi ++; sz <<= 1; } ll inv_sz = inv(sz); static ll tmp[maxn]; std::copy(B, B + mod, tmp); std::fill(tmp + mod, tmp + sz, 0LL); NTT(tmp, bi); NTT(BB, bi); for ( int i = 0; i < sz; i ++) { tmp[i] = (tmp[i] * BB[i]) % ha; tmp[i] = (tmp[i] * (ha - 1LL)) % ha; tmp[i] = (tmp[i] + 2LL) % ha; tmp[i] = (tmp[i] * BB[i]) % ha; } NTT(tmp, bi, true ); for ( int i = 0; i < sz; i ++) { tmp[i] = (tmp[i] * inv_sz) % ha; } std::copy(tmp, tmp + mod, BB); std::fill(BB + mod, BB + sz, 0LL); } } int main() { static ll fac[maxn], A[maxn], B[maxn], BB[maxn]; int n; scanf ( "%d" , &n); int bi = 0, sz = 1; while (sz <= n + 1) { bi ++; sz <<= 1; } fac[0] = 1LL; for ( int i = 1; i <= n; i ++) { fac[i] = (fac[i - 1] * (ll(i))) % ha; } B[0] = 1LL; for ( int i = 1; i <= n; i ++) { B[i] = pow_mod(2LL, (ll(i)) * (ll(i - 1)) / 2LL); B[i] = (B[i] * inv(fac[i])) % ha; } for ( int i = 1; i <= n; i ++) { A[i] = pow_mod(2LL, (ll(i)) * (ll(i - 1)) / 2LL); A[i] = (A[i] * inv(fac[i - 1])) % ha; } poly_inv(n + 1, B, BB); poly_mul(A, BB, bi, A); printf ( "%lld\n" , (A[n] * fac[n - 1]) % ha); return 0; } |