不妨设\(A(x)\)表示答案。
对于一个点，考虑它的三个子节点，直接卷起来是\(A(x)^3\)，但是这样肯定会计重，因为我们要的是无序的子节点。
那么用burnside引理，枚举一个排列，一个环的选择要相同，如果环的大小是y，则对应\(A(x^y)\)。
最后可以得到：
\(A(x)=x{A(x)^3+3A(x^2)A(x)+2A(x^3)\over 6}+1\)
分治NTT可以解这个方程，不过因为有3次的，比较复杂，考虑用牛顿迭代：
\(F(A(x))=x{A(x)^3+3A(x^2)A(x)+2A(x^3)\over 6}+1-A(x)=0\)
\(new A(x)=A(x)-{F(A(x))\over F(A(x))'}\)
问题在于\(F(A(x))\)中有\(A(x^2)、A(x^3)、x\)这样的项，如何求导。
注意牛顿迭代的倍增中，已经求出了\(mod~x^{n/2}\)的答案，那么\(A(x^2)、A(x^3)、x\)是已知的，可以视作常数。
所以\(F(A(x))'=x{3A(x)^2 + 3A(x^2)\over 6} - 1\)。
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i < B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;
const int mo = 998244353;
ll ksm(ll x, ll y) {
ll s = 1;
for(; y; y /= 2, x = x * x % mo)
if(y & 1) s = s * x % mo;
return s;
}
typedef vector<ll> V;
#define pb push_back
#define si size()
#define re resize
namespace ntt {
const int nm = 1 << 18;
int r[nm]; ll w[nm], a[nm];
void build() {
for(int n = 1; n < nm; n *= 2) {
ll v = ksm(3, (mo - 1) / 2 / n) ;
w[n] = 1;
ff(i, 1, n) w[n + i] = w[n + i - 1] * v % mo;
}
}
void dft(V &c, int f) {
int n = c.si;
ff(i, 0, n) a[i] = c[i];
ff(i, 0, n) {
r[i] = r[i / 2] / 2 + (i & 1) * n / 2;
if(i < r[i]) swap(a[i], a[r[i]]);
} ll b;
for(int i = 1; i < n; i *= 2) for(int j = 0; j < n; j += 2 * i) ff(k, 0, i)
b = a[i + j + k] * w[i + k], a[i + j + k] = (a[j + k] - b) % mo, a[j + k] = (a[j + k] + b) % mo;
if(f == -1) {
reverse(a + 1, a + n);
b = ksm(n, mo - 2);
ff(i, 0, n) a[i] = (a[i] + mo) * b % mo;
}
ff(i, 0, n) c[i] = a[i];
}
}
using ntt :: dft;
V operator *(V a, V b) {
int n0 = a.si + b.si - 1, n = 1;
while(n < n0) n *= 2;
a.re(n); b.re(n);
dft(a, 1); dft(b, 1);
ff(i, 0, n) a[i] = a[i] * b[i] % mo;
dft(a, -1);
a.re(n0); return a;
}
V operator +(V a, V b) {
a.re(max(a.si, b.si));
ff(i, 0, b.si) a[i] = (a[i] + b[i]) % mo;
return a;
}
V operator -(V a, V b) {
a.re(max(a.si, b.si));
ff(i, 0, b.si) a[i] = (a[i] - b[i] + mo) % mo;
return a;
}
V operator * (V a, int b) {
ff(i, 0, a.si) a[i] = a[i] * b % mo;
return a;
}
V qni(V a) {
V b; b.resize(1); b[0] = ksm(a[0], mo - 2);
for(int n = 2; n < a.si * 2; n *= 2) {
V c = a; c.re(n); c.re(2 * n);
V d = b; b.re(2 * n); d.re(n);
dft(b, 1); dft(c, 1);
ff(i, 0, b.si) b[i] = b[i] * b[i] % mo * c[i] % mo;
dft(b, -1); b.re(n);
ff(i, 0, n) b[i] = (2 * d[i] - b[i] + mo) % mo;
}
b.re(a.si); return b;
}
V yy(V a) {
a.insert(a.begin(), 0);
return a;
}
V stay(V a, int y, int n) {
V b; b.resize(n);
int mx = a.si - 1; mx = min(mx, (n - 1) / y);
fo(i, 0, mx) b[i * y] = a[i];
return b;
}
const ll ni6 = ksm(6, mo - 2);
V newton(int n0) {
V b; b.re(1); b[0] = 1;
for(int n = 2; n <= n0 * 2; n *= 2) {
V c = b * b * b + stay(b, 2, n) * b * 3 + stay(b, 3, n) * 2;
c = yy(c * ni6); c.re(n);
c[0] ++; c = c - b;
V d = b * b * 3 + stay(b, 2, n) * 3;
d = yy(d * ni6); d[0] --; d.re(n);
b = b - c * qni(d); b.re(n);
}
return b;
}
int n;
int main() {
ntt :: build();
scanf("%d", &n);
V a = newton(n);
pp("%lld\n", a[n]);
}