Asusetic
0%

烷基计数

发表于 分类于 Valine：
本文字数： 1.8k 阅读时长 ≈ 2 分钟

Polya的奇怪应用

Update：真的polya定理

烷基计数

参考：神仙的化学计数

给出一个比较通用的计数方法

定义Z(G)Z(G)为有限置换群GG的循环指数，Z(G)=1∣G∣∑g∈G∏rsrcr(g)Z(G)=\frac{1}{|G|}\sum\limits_{g\in G}\prod\limits_rs_r^{c_r(g)}，这是一个关于sis_i的多项式

cr(g)c_r(g)是gg的循环分解中有rr个元素的循环，也就是将gg循环分解表示加起来

比如Z(S3)=16(2s3+3s2s1+s13)Z(S_3)=\frac{1}{6}(2s_3+3s_2s_1+s_1^3) (如果考虑手性就是群A3\mathfrak{A}_3)

而如果我们把sis_i用多项式Z[xr]\mathbb{Z}[x^r]代替就是Z(S3,F(x))=16(2F(x3)+3F(x2)F(x)+F3(x))Z(S_3,F(x))=\frac{1}{6}(2F(x^3)+3F(x^2)F(x)+F^3(x)) （这其实就是真的Polya定理

考虑烷基计数

拆去第一个碳原子，剩下的取代基可以自由分配，所以其生成函数F(x)=1+xZ(s3,F(X))=1+x6[2F(x3)+3F(x2)F(x)+F3(x)]F(x)=1+xZ(s_3,F(X))=1+\frac{x}{6}[2F(x^3)+3F(x^2)F(x)+F^3(x)]

可以分治NTT或者牛顿迭代

选择不用脑子的牛顿迭代

G(F(X))=1+x6[2F(x3)+3F(x2)F(x)+F3(x)]−F(x)G(F(X))=1+\frac{x}{6}[2F(x^3)+3F(x^2)F(x)+F^3(x)]-F(x)

设已经求出H(x)H(x)满足G(H(x))≡0(modxn2)G(H(x))\equiv0\pmod{ x^{\frac{n}{2}}}

那么F(x)=H(x)−G(H(x))G′(H(x))F(x)=H(x)-\frac{G(H(x))}{G'(H(x))}

G′(H(x))G'(H(x))中因为已知可以当作常数G′(H(x))=x3H(x)2+3H(x2)6−1G'(H(x))=x\frac{3H(x)^2+3H(x^2)}{6}-1

带入F(X)=H(x)−6−x(H(x)3+3H(x2)H(x)+2H(x3))−6H(x)3x(H(x)2+H(x2))−6F(X)=H(x)-\frac{6-x(H(x)^3+3H(x^2)H(x)+2H(x^3))-6H(x)}{3x(H(x)^2+H(x^2))-6}

实际上这个方法像烯烃/烷烃计数都通用的，但是我懒得再写了

Tips:烷烃的话答案是

G(x)=x[124F4(x)+14F2(x)F(x2)+18F2(x2)+13F(x)F(x3)+14F(x4)]−12F2(x)+12F(x2)+F(x)G(x)=x[\frac{1}{24}F^4(x)+\frac 14 F^2(x)F(x^2)+\frac{1}{8}F^2(x^2)+\frac{1}{3}F(x)F(x^3)+\frac 14 F(x^4)]-\frac 12 F^2(x)+\frac 12 F(x^2)+F(x)

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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
#include <algorithm>
// #include <dbg>
#include <iostream>
#include <vector>
using namespace std;
namespace Poly {
typedef long long ll;
const ll mod = 998244353;
const ll G = 3, invG = 332748118;
typedef vector<ll> poly;
const int N = 4e6 + 15;
int rev[N];
ll qpow(ll a, int b) {
    ll res = 1;
    while (b) {
        if (b & 1) {
            res = res * a % mod;
        }
        a = a * a % mod;
        b /= 2;
    }
    return res;
}
void NTT(poly &a, int len, int typ) {
    a.resize(len);
    int lim = a.size();
    for (int i = 0; i < lim; i++) {
        if (i < rev[i]) {
            swap(a[i], a[rev[i]]);
        }
    }
    for (int mid = 1; mid < lim; mid *= 2) {
        ll wn = qpow(typ == 1 ? G : invG, (mod - 1) / (mid * 2));
        for (int j = 0; j < lim; j += mid * 2) {
            ll w = 1;
            for (int k = 0; k < mid; k++) {
                ll x = a[j + k], y = w * a[j + k + mid] % mod;
                a[j + k]       = (x + y) % mod;
                a[j + k + mid] = (x - y + mod) % mod;
                w              = w * wn % mod;
            }
        }
    }
    if (typ == 1) {
        return;
    }
    int invdeg = qpow(lim, mod - 2);
    for (int i = 0; i < lim; i++) {
        a[i] = a[i] * invdeg % mod;
    }
}
ll W[N];
void prework(int n) {
    for (int i = 1; i < n; i <<= 1) {
        W[i]  = 1;
        int w = qpow(3, (mod - 1) / (i * 2));

        for (int j = 1; j < i; j++)
            W[i + j] = (ll)W[i + j - 1] * w % mod;
    }
}

int getLen(int deg) {
    int len = 1;
    while (len < deg) {
        len *= 2;
    }
    return len;
}
void getRev(int deg) {
    for (int i = 0; i < deg; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (deg >> 1) : 0);
    }
}
int PrepareNTT(int deg) {
    int L = getLen(deg);
    getRev(L);
    return L;
}
poly _mul(const poly &a, const poly &b) {
    int len = a.size();
    poly res(len);
    for (int i = 0; i < len; i++) {
        res[i] = a[i] * b[i] % mod;
    }
    return res;
}
poly operator-(const poly &a) {
    int len = a.size();
    poly res(len);
    for (int i = 0; i < len; i++) {
        res[i] = (mod - a[i]) % mod;
    }
    return res;
}
poly operator+(const poly &a, const poly &b) {
    int len = a.size();
    poly res(len);
    for (int i = 0; i < len; i++) {
        res[i] = a[i] + b[i] % mod;
    }
    return res;
}
poly operator-(const poly &a, const poly &b) {
    return a + -b;
}
poly inv(poly a) {
    poly res, inv_a;
    res.resize(1);
    res[0]  = qpow(a[0], mod - 2);
    int lim = getLen(a.size());
    for (int len = 2; len <= lim; len *= 2) {
        int n = PrepareNTT(len + len);
        inv_a = a;
        inv_a.resize(n);
        res.resize(n);
        for (int i = len; i < n; i++) {
            inv_a[i] = 0;
        }
        NTT(inv_a, n, 1);
        NTT(res, n, 1);
        for (int i = 0; i < n; i++) {
            res[i] = 1ll * (2ll - 1ll * inv_a[i] * res[i] % mod + mod) % mod *
                     res[i] % mod;
        }
        NTT(res, n, -1);
        for (int i = len; i < n; i++) {
            res[i] = 0;
        }
    }
    res.resize(a.size());
    return res;
}
} // namespace Poly
using namespace Poly;
ll solve(int n) {
    poly c(1, 1), hx2, hx3, p, q;
    for (int len = 1; len < n; len *= 2) {
        hx2 = hx3 = poly(len * 2);
        for (int j = 0; j < len; j++)
            hx2[j << 1] = c[j];
        for (int j = 0; j * 3 < 2 * len; j++)
            hx3[j * 3] = c[j];
        PrepareNTT(len * 4);
        NTT(c, len * 4, 1);
        NTT(hx2, len * 4, 1);
        NTT(hx3, len * 4, 1);
        p = q = poly(len * 4);
        for (int j = 0; j < (len * 4); j++) {
            int w = j < (len * 2) ? W[len * 2 | j] : mod - W[len * 2 | j];
            q[j] =
                (3ll * w * ((ll(c[j]) * c[j] % mod + hx2[j]) % mod) + mod - 6) %
                mod;
            p[j] = (w *
                        ((ll(c[j]) * c[j] % mod * c[j] +
                          3ll * hx2[j] * c[j] % mod + 2ll * hx3[j]) %
                         mod) %
                        mod +
                    (mod - 6ll) * c[j] % mod + 6ll) %
                   mod;
        }

        NTT(p, len * 4, -1);
        NTT(q, len * 4, -1);
        q.resize(len * 2);
        q = inv(q);
        p.resize(len * 2);
        NTT(p, len * 4, 1);
        NTT(q, len * 4, 1);
        for (int j = 0; j < (len * 4); j++) {
            c[j] = (c[j] + (mod - p[j]) * (ll)q[j]) % mod;
        }
        NTT(c, len * 4, -1);
        c.resize(len * 2);
    }
    return c[n - 1];
}
int n;
poly res;
int main() {
    prework(2e6 + 10);
    cin >> n;
    cout << solve(n + 1);
    return 0;
}