Prufer编码 & [bzoj 1005] [HNOI2008]明明的烦恼：Prufer编码，组合数学，高精度

题意：明明是个奇怪的孩子，想知道n个有标号结点（0<n≤1000）、一些结点的度数已经给定的生成树有多少棵，还不让你取模。

以前听说n个点的完全图的生成树有$n^{n-2}$棵，从Matrix67前辈的博客中学到用prufer编码证明这个结论，并且看到一个关于“每个结点度数已知的生成树”棵数的公式。

读了题，联想到那些。然而我已经忘得差不多了。没有立即重新学习，在记忆中搜寻，运用一下排列组合知识，发挥一下创造力，终于脑补出来了……

不考虑一个结点的生成树。两个结点的生成树只有1棵，定义它的prufer编码为空。称度数为1的结点为叶子。n≥3时，我们这样求prufer编码（简洁起见，递归地描述）：

找叶子中标号最小的那一个，把和它相邻的点记下，删掉这片叶子，求剩下这棵树的prufer编码。n=2时停止。

我们得到一个长度为(n-2)的序列。每棵树对应的序列是唯一的。那么，每个序列对应的树是否唯一？

叶子从不出现在prufer序列中。设点u和叶子v相邻，当u为叶子时，图中必然只剩u、v两点，根据定义，prufer编码为空。

不出现在prufer序列中的点一定是叶子。对于终被删除且从未在序列中出现过的点v，它被删除的时候一定是叶子；由于它没在序列中出现，它不曾和已被删除的点相邻，也就是说，它成为叶子并不是因为相邻的点被删除，而是与生俱来。对于n=2时的幸存者v，这时v也是叶子，论证同理。

于是，我们得到一个解码过程。

取没在prufer序列中出现过的标号最小的点，把它和序列中的第一个点相连，删掉这个点，对从第二项起的序列解码。序列为空时，将剩下的两个点相连。

最终一定剩下两个点，因为删除了(n-2)个点。由于不出现在序列中的点就是叶子，我们只是把prufer编码的过程重做了一遍，不同之处在于，编码是到图中看点的度数，解码是在序列中看哪些点没出现过。

树对应唯一的编码，每段编码有且仅有一棵树与之对应（如果两棵树T、T’编码后都得到x，x解码既得到T又得到T’，推出T=T’），因此，这是一个一一映射，n个点的完全图的生成树有$n^{n-2}$棵。

如果v在序列中出现(d-1)次，那么v的度数为d。每个点最终都要成为叶子，成为叶子之后它就不会再在序列中出现。除了与作为叶子的v相邻的点u，其他原先与v相邻的点要么不存在（d-1=0，d=1），要么被删去，每次删除时v都出现一次，总共(d-1)次，加上u，v的度数为d。

回到本题。设$v_1, v_2, \cdots, v_k(0 \le k \le n)$的度数被限制为$d_1, d_2, \cdots, d_k$，先安排$v_i(1 \le i \le k)$在序列中的位置（可重集的排列问题），在剩下的空位中不限次数地选择剩余的(n-k)个结点，即得答案：

$$answer = \frac {(n-2)!(n-k)^r} {r!(d_1-1)!(d_2-1)!\cdots(d_k-1)!}, r = n-2-\sum_{i=1}^k(d_i-1)$$

俗话说得好，计数问题不取模就是耍流氓。本题中明明在耍流氓。为了避免高精度除法，我们按质因子保存，最后再来个高精度乘法。

然后开始犯蠢：
1. 我竟然认为n的质因子最大为$\sqrt n$。
2. 没特判n=1。
3. 我的高精度没保存每个数实际使用的长度，每次334*334地暴力乘，导致TLE。

针对问题3进行修正后不超时了，然而240ms还是比别人慢好多……原因应该是我预处理了1~n阶乘的质因数分解，这样可以递推，要用的时候直接取。事实上只会用到不超过(k+2)个；每次虽然得把不大于x的数都分解一遍，但有可能数据中的x较小。

关于阶乘的质因数分解，hzwer的博客中提到一种方法（并注明转自怡红公子），可以参考这篇文章：阶乘质因数分解。

#include <cstdio>
#include <cstring>
const int MAX_N = 1000, M = 168;
using namespace std;
typedef long long ll;
int n, prime[M] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,
    199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,
    457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,
    733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997};

struct Fac {
    int x[M];
    int& operator[](int i)
    {
        return x[i];
    }
    Fac operator+=(Fac& b)
    {
        for (int i = 0; i < M; ++i)
            x[i] += b[i];
        return *this;
    }
    Fac operator-=(Fac& b)
    {
        for (int i = 0; i < M; ++i)
            x[i] -= b[i];
        return *this;
    }
} f[MAX_N+1], g;

struct Big {
    static const int w = 334, lg = 9, base = 1e9;
    int x[w], l;
    Big(int a=0)
    {
        *this = a;
    }
    Big operator=(int a)
    {
        l = 1;
        x[0] = a;
    }
    Big operator=(const Big& b)
    {
        memcpy(x, b.x, sizeof(int)*(l = b.l));
    }
    Big operator*=(const Big& b)
    {
        Big c;
        c.l = w<l+b.l ? w : l+b.l;
        memset(c.x, 0, sizeof(int)*c.l);
        for (int i = 0; i < l; ++i)
        for (int j = 0, f = 0; i+j < c.l; ++j) {
                ll t = c.x[i+j] + (ll)x[i]*(j < b.l ? b.x[j] : 0) + f;
                c.x[i+j] = t % base;
                f = t / base;
        }
    if (c.l > 1 && !c.x[c.l-1])
        --c.l;
        return *this = c;
    }
    void out(char c='\n')
    {
        int p = l-1;
        while (p && !x[p])
            --p;
        printf("%d", x[p]);
        while (p)
            printf("%0*d", lg, x[--p]);
        putchar(c);
    }
    friend Big pow(Big x, int n)
    {
        Big y(1);
        while (n) {
            if (n & 1)
                y *= x;
            x *= x;
            n >>= 1;
        }
        return y;
    }
};

void init()
{
    for (int i = 2; i <= n; ++i)
        for (int j = 0; j < M; ++j)
            if (i % prime[j] == 0) {
                f[i][j] = 1;
                f[i] += f[i/prime[j]];
                break;
            }
    for (int i = 3; i <= n; ++i)
        f[i] += f[i-1];
}

int main()
{
    scanf("%d", &n);
    if (n == 1) {
        int d;
        scanf("%d", &d);
        printf("%d\n", d <= 0);
        return 0;
    }
    init();
    g = f[n-2];
    int sum = 0, cnt = 0;
    for (int i = 0; i < n; ++i) {
        int d;
        scanf("%d", &d);
        if (!d) {
            puts("0");
            return 0;
        }
        if (--d >= 0) {
            g -= f[d];
            sum += d;
            ++cnt;
        }
    }
    int b = n-2-sum;
    if (b < 0) {
        puts("0");
        return 0;
    } 
    g -= f[b];

    Big ans(pow(Big(n-cnt), b));
    for (int i = 0; i < M; ++i)
        if (g[i])
            ans *= pow(Big(prime[i]), g[i]);

    ans.out();

    return 0;
}