这道题的主要目的是普及无标号无根树的计数方法，故不讨论非多项式时间做法。其实 $n$ 完全可以出得更大，但是没必要。
这道题改自 ProjectEuler#677。
思路要点
首先考虑无标号有根树，我们只需维护 dp 数组 $R(n), B(n), Y(n)$ 即可。通过更新 dp 数组 $Q(d, n)$ 表示此时不在意节点的根的颜色，且根的度数为 $d$ 即可维护红根和蓝根，$W(d, n)$ 用于维护孩子没有黄根树的，这可以帮助计算黄根。每次 update 会消耗 $\Theta(n)$ 的时间，比如说已经计算完了 $R(n)$ 的值，那么它作为子树有可能出现了 $k$ 次，则选出 $\displaystyle \binom{R(n) + k - 1}{k - 1}$ 种方案，令 $Q(d,N)$ 加上 $\binom{R(n) + k - 1}{k - 1} Q(d - k, N - kn)$。
这一部分的复杂度为 $\Theta(n^2)$。
计算无标号无根树的要点是观察树的一个特殊点：一颗树要么存在唯一一个重心使得任何一个子树的大小均 $< \frac{n}2$，要么存在一个“重边”使得它将树恰好分为大小为 $\frac n2$ 的两部分。
对于前者部分的计数，我们只需将 dp 数组计算到 $\lceil \frac n2 \rceil - 1$ 时，将 $Q,R$ 数组取出。对于后者，若 $n$ 是偶数，当计算完 $\frac n2$ 的 dp 值时，将 $R, B, Y$ 用于计算即可。
如何做到更快
考虑计算 $R,B,Y$ 的生成函数 $G_R(x), G_B(x), G_Y(x)$。
我们记 $S = G_R + G_B + G_Y$，我们使用 Pólya 计数定理表示 $3$ 个孩子的树，通过枚举置换群：
$$
\frac{S(x)^3 + 3S(x^2)S(x) + 2S(x^3)}{3!}
$$
读者可自己尝试推导 $4$ 和 $2$ 个孩子的树。
由此我们可以通过 $G_R, G_B, G_Y$ 自身带入 $x, x^2, x^3$ 这些的多项式来表示。接下来考虑如何求解。
我们让 $G_R, G_B, G_Y$ 放在一起进行迭代，我们考虑倍增，由于已经计算出了前 $\frac n2$ 项，所有带入 $x^2, x^3$ 的部分就已经确定了，我们可以认为是常数。接下来就是一个关于 $G_R, G_B, G_Y$ 的三元三次方程，这依然可以通过牛顿迭代解决，也就是这等价于根据 $G_R, G_B, G_Y$ 的定义式得到了 $3$ 个它们之间的关系式 $E_R/E_B/E_Y(G_R, G_B, G_Y) = 0$，将 $E_R/E_B/E_Y(u, v, w)$ 在 $G_R, G_B, G_Y$ 处泰勒展开就是一次方程，求解一个三元一次方程即可。
本算法的复杂度为 $\Theta(n\log n)$，常数显而易见的大，欢迎有兴趣的同学来实现。
事实上本题是 UOJ #450. 【集训队作业2018】复读机 的一个加强版。
算法一
$n\le 5\times 10^4$，我可以做多项式快速幂！这个多项式就是 $(\sum_{j\ge 0} \frac{x^{jd}}{(jd)!})^k$。这个模数不太友好，你可能需要 MTT……常数应该很大。
时间复杂度 $\Theta(n\log n) \sim \Theta(n\log n \log k)$，前者是进行多项式 $\ln , \exp$ 达到的复杂度，但是常数可能很大。
预计得分 $20\%$。
算法二
我们考察 $(\sum_{j\ge 0} \frac{x^{jd}}{(jd)!})^k$ 可以如何表示。
事实上与周期函数有关的，我们可以自然想到类似 DFT 的形式，也就是考虑到 $\sum_{j=0}^{d - 1} \omega_d^{cj} = d[d | c]$，因此我们可以知道上面的那个多项式实际上是
$$
\frac1d \sum_{j=0}^{d-1} \exp \omega_d^j x
$$
我们直接枚举拆括号的过程中计算了几个 $\exp \omega_d^0 x$，几个 $\exp \omega_d^1 x$……然后这一部分的贡献就是 $(c_0 + c_1 \omega_d^1 + c_2 \omega_d^2 + \dots)^n$。最后总和除以 $d^k$。
时间复杂度 $\Theta(k^{d-1}\log n)$，预计得分 $50\%$。
算法三
主要思想还是围绕优化计算下式：
$$
[x^n]\left( \frac1d \sum_{j=0}^{d-1} \exp \omega_d^j x \right)^k
$$
对于 $d=4$ 的情况，我们实际上是只需要统计每个 $\exp (a + b\mathrm{i})x$ 出现多少次。通过推导一些式子，或者直接将复平面旋转 $45^\circ$ 我们可以得到，方案数是 $\displaystyle\binom{k}{\frac{k + |a+b|}2} \binom{k}{\frac{k + |a-b|}2}$。因此答案可以在 $\Theta(k^2 \log n)$ 内计算出来。
预计得分 $70\%$。
算法四
对于 $d=6$ 的情况，我们注意到 $\omega_6$ 的代数关系有 $\omega - 1 = \omega^2$，因此所有的根都可以用 $1, \omega$ 表示，即给每个数赋予了一个离散的二维坐标
$$
\begin{array}{clclc}
\omega^0 & = & 1 & & \\
\omega^1 & = & & & \omega\\
\omega^2 & = & -1 & + & \omega\\
\omega^3 & = & -1 & & \\
\omega^4 & = & & &-\omega\\
\omega^5 & = & 1 & + &-\omega
\end{array}
$$
因此我们只需要做一个二维的快速幂，倍增 FFT 可以做到 $\Theta(k^2 \log k)$。但是常数可能较大。
预计得分 $80\% \sim 100\%$。
算法五
考虑计算一个二元母函数的高阶幂，$G(x, y) = F(x, y)^k$，可以得到 $F \frac{\partial G}{\partial_x} = kG \frac{\partial F}{\partial_x}$，由此可列得递推式子解出所有系数。计算高阶幂的复杂度是 $\Theta(k^2)$ 的。复杂度 $\Theta(k^2 \log n)$ 且常数较小。
事实上这也是可以直接用来做 $k=4$ 的情况的。
预计得分 $100\%$。
further
对于更大的 $d$，我们期望能得到一个怎样的做法？考虑 $d$ 次单位根的代数关系，我们希望能够基于一个它们的有理系数线性组合的基。由此则能将维度缩到基的维度，在其上进行快速幂。
我将说明维度可以达到 $\varphi(d)$，而数学迷告诉我说通过一些关于域的理论可以证明这也是下界，等我学了之后补一下证明。
考虑分圆多项式 $\Phi_d(x) = \prod_{0\le k < d, \gcd(d, k) = 1} (x - \omega_d^k)$，它的次数显然是 $\varphi(d)$。
显然分圆多项式也可以用容斥原理重写，即
$$
\Phi_d(x) = \prod_{k|d} (x^{d/k}-1)^{\mu(k)}
$$
显然该多项式的系数都是整数。由于 $\Phi_d(\omega_d) = 0$，不难得到 $\omega_d^k = \left.x^k \bmod \Phi_d\right \vert_{x=\omega}$。因此这个多项式的取模自然而然地导出了每个单位根到 $1, \omega_d, \dots, \omega_d^{\varphi(d) - 1}$ 的线性组合。
算法一
显然所选的圆上要么有两个点要么有至少三个。因此分别枚举，枚举两个点作为圆的直径，此外枚举三个点之后可以确定唯一一个外接圆。
时间复杂度 $\Theta(n^4)$，预计得分 $10\%$。
算法二
$m=n$ 的情况就是经典的最小圆覆盖问题。
存在期望时间复杂度 $\Theta(n)$ 的算法，预计得分 $20\%$。
算法三
我们考虑二分答案。那么接下来枚举哪个点在圆上，其它每个点我们可以计算出它对于当前圆在哪个夹角区间内，进行一遍前缀和即可查找出是否存在一个夹角使得该圆包含 $m$ 个点。
时间复杂度 $\Theta(n^2\log n \log \frac{x}{\epsilon})$，预计得分 $60\%$。
算法四
我们对算法三进行进一步的观察。事实上我们可以看成这样一个问题，对于第 $i$ 个点，客观上存在一个 $r_i$ 即最小的圆的半径使得圆内有 $m$ 个点（我们认为无解即 $+\infty$），我们在二分的过程中，总共会预计调用 $\Theta(n\log \frac{x}{\epsilon})$ 次检查。这实际上是有所浪费的，在同样的检查次数内，实际上我们还可以算出每个 $i$ 对应的 $r_i$，换句话说我们多获取了不必要的信息，因为我们只关心最小值。
接下来就是有趣的地方了：我们先将序列随机打乱，并把当前最小值设为 $a = +\infty$。接下来我们将点一个个考虑，我们实际上可以一次检查出当前点的答案是否可以 $< a$。如果并不小于，我们就可以略过它，否则对该点的答案进行二分。这样显然是进行 $n + L\log \frac{x}{\epsilon}$ 次检查，其中 $L$ 是 $r_i$ 序列的单调栈长度。
显然最坏情况是 $r$ 互不相同时，单调栈长度的期望最长，这等价于一个排列的期望长度。而这一期望长度等于 $H_n = \sum_{k=1}^n \frac 1k = \Theta(\log n)$。
因此本算法的期望复杂度为 $\Theta( (n + \log n \log \frac{x}{\epsilon}) n\log n)$，预计得分 $100\%$。
