洛谷P4727-HNOI2009-图的同构计数

题解写得很好：https://www.luogu.com.cn/blog/pythoner713/solution-p4727

但是Part 5不是很清晰。这里用自己的语言重述Part 5。

Part 4的结论是，对于一个置换\(g\)，如果其循环表示的循环的长度分别为\(b_1, b_2, \cdots, b_K\)，那么在这个置换下，边的等价类的数目

\[k = \sum_{i=1}^K \lfloor \frac{b_i}{2} \rfloor + \sum_{i=1}^K \sum_{j=1}^{i-1} gcd(b_i, b_j) \]

\(g\)的不动点满足每个等价类里的边都要染成同一个颜色，所以由Burnside引理，不同构的染色方案数为

\[\frac{1}{|G|} \sum_{g\in G} 2^k \]

可见，每个置换对答案的贡献仅仅取决于其循环的长度的多重集\([b_1, b_2, \cdots, b_K]\)。我们可以枚举每个这样的循环长度的集合，算出对应的置换的个数，这样就可以把这些对应的置换对答案的贡献一次性算出来了。

对于这样一个循环长度的多重集：\([b_1, b_2, \cdots, b_K]\)，先将每个循环看作是不同的，将每个点分配到一个循环里，分配方案有\(\frac{n!}{\prod_{i=1}^K b_i!}\)种。然后每个循环内部要将这些点组织成一个循环，相当于一个圆排列，有\((b_i - 1)!\)种方式。此外，相同的长度的循环本质上是相同的，假设有\(C\)种不同长度：\(\{B_1, B_2, \cdots, B_C\}\)，长度为\(B_i\)的循环有\(c_i\)个，那么答案还要除以\(\prod_{i=1}^C c_i!\)。循环长度的多重集\([b_1, b_2, \cdots, b_K]\)对应的置换的个数为

\[\frac{n!}{\prod_{i=1}^K b_i!} \prod_{i=1}^K(b_i - 1)! \frac{1}{\prod_{i=1}^C c_i!} = \frac{n!}{\prod_{i=1}^K b_i \prod_{i=1}^C c_i!} \]

将枚举的循环长度的多重集记为\(b\)，其对应的边的等价类个数为\(k\)，那么答案为

\[\frac{1}{|G|} \sum_{b} \frac{n!}{\prod_{i=1}^K b_i \prod_{i=1}^C c_i!} 2^k \]

由于置换个数\(|G| = n!\)，所以最终的式子为

\[\sum_{b} \frac{2^k}{\prod_{i=1}^K b_i \prod_{i=1}^C c_i!} \]

枚举循环长度的多重集实际上是枚举\(n\)的整数拆分。60的整数拆分有966467个。每个循环长度的多重集对答案的贡献可以在\(O(K^2)\)求出。所以总的时间大概是\(1e8\)的样子。

rust代码：

use std::io;

const P: u32 = 997;

// Faster?
fn mpow(mut base: u32, mut ex: u32, p: u32) -> u32 {
    let mut ans = 1u32;
    while ex != 0 {
        if ex & 1u32 != 0 {
            ans = (ans as u64 * base as u64 % p as u64) as u32;
        }
        base = (base as u64 * base as u64 % p as u64) as u32;
        ex >>= 1;
    }
    return ans;
}

fn get_invs(invs: &mut [u32], p: u32, n: u32) {
    invs[1] = 1;
    for i in 2..(n + 1) {
        invs[i as usize] = (((p - p / i) as u64) * (invs[(p % i) as usize] as u64) % p as u64) as u32;
    }
}

fn gcd(mut a: u32, mut b: u32) -> u32 {
    while b != 0 {
        let r = a % b;
        a = b;
        b = r;
    }
    return a;
}

fn work(b: &Vec<u32>, invs: &[u32; P as usize]) -> u32 {
    let mut k = 0;
    for bi in b {
        k += bi / 2;
    }
    for i in 0..b.len() {
        for j in 0..i {
            k += gcd(b[i], b[j]);
        }
    }
    let mut ans = mpow(2, k, P);
    for bi in b {
        ans = ans * invs[*bi as usize] % P;
    }
    let mut last = 0;
    let mut ci_fact = 1;
    let mut ci = 0;
    for bi in b {
        let bi = *bi;
        if bi == last {
            ci += 1;
            ci_fact = ci_fact * ci % P;
        } else {
            ans = ans * invs[ci_fact as usize] % P;
            last = bi;
            ci = 1;
            ci_fact = 1;
        }
    }
    ans = ans * invs[ci_fact as usize] % P;
    return ans;
}
fn dfs(n: u32, mn: u32, b: &mut Vec<u32>, inv: &[u32; P as usize]) -> u32 {
    if n == 0 {
        return work(b, inv);
    } else if n < mn {
        return 0;
    }
    let mut ans = 0;
    for bi in mn..(n+1) {
        b.push(bi);
        ans = (ans + dfs(n - bi, bi, b, inv)) % P;
        b.pop();
    }
    ans
}

fn main() {
    let mut input = String::new();
    io::stdin().read_line(&mut input).unwrap();
    let n = input.trim().parse().unwrap();
    let mut b = Vec::new();
    let mut invs = [0; P as usize];
    get_invs(&mut invs, P, P - 1);
    println!("{}", dfs(n, 1, &mut b, &invs));
}