【題目連結】

• 點選開啟連結

【思路要點】

• 考慮用 $Burnsid$
  e Burnside 引理計數。
• 不妨令 $N≤M$
  ，列舉 $N$ 的整數拆分，在第二維上 $dp$
  即可。
• 具體來說，滿足將 $N$ 拆分為 $N=\sum a_ic_i\ (a_i>a_{i-1})$ 的置換數為 $\frac{N!}{\prod a_i^{c_i}*c_i!}$ ，一對環長為 $(a,b)$ 的置換環對不動點個數產生的貢獻為 $*2^{gcd(a,b)}$ ，列舉 $N$ 的拆分後，可以通過簡單 $dp$ 計算所需貢獻。
• 令 $N≤M$ ，時間複雜度 $O(f(N)*M^2LogM)$ ，其中 $f(N)$ 表示 $N$ 的整數拆分數。

【程式碼】

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 1005;
const int P = 1e9 + 7;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
template <typename T> void chkmax(T &x, T y) {x = max(x, y); }
template <typename T> void chkmin(T &x, T y) {x = min(x, y); } 
template <typename T> void read(T &x) {
	x = 0; int f = 1;
	char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;
	for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';
	x *= f;
}
template <typename T> void write(T x) {
	if (x < 0) x = -x, putchar('-');
	if (x > 9) write(x / 10);
	putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) {
	write(x);
	puts("");
}
int fac[MAXN], inv[MAXN], two[MAXN * MAXN], g[MAXN][MAXN];
int power(int x, int y) {
	if (y == 0) return 1;
	int tmp = power(x, y / 2);
	if (y % 2 == 0) return 1ll * tmp * tmp % P;
	else return 1ll * tmp * tmp % P * x % P;
}
int n, m, cnt, ans, a[MAXN], dp[MAXN][MAXN];
void update(int &x, int y) {
	x += y;
	if (x >= P) x -= P;
}
int getdp() {
	for (int i = 0; i <= m; i++)
	for (int j = 0; j <= m; j++)
		dp[i][j] = 0;
	dp[m][m] = 1;
	for (int i = m; i >= 0; i--)
	for (int j = m; j >= 1; j--) {
		if (dp[i][j] == 0) continue;
		int tmp = 0;
		for (int k = 1; k <= cnt; k++)
			tmp += g[j][a[k]];
		int val = dp[i][j];
		for (int k = 0; j * k <= i; k++) {
			update(dp[i - j * k][j - 1], 1ll * val * inv[k] % P);
			val = 1ll * val * inv[j] % P * fac[j - 1] % P * two[tmp] % P;
		}
	}
	return 1ll * dp[0][0] * fac[m] % P;
}
void check(int depth) {
	cnt = depth;
	int now = fac[n];
	for (int i = 1; i <= cnt; i++) {
		now = 1ll * now * inv[a[i]] % P;
		now = 1ll * now * fac[a[i] - 1] % P;
		if (a[i] != a[i - 1]) {
			int j = i;
			while (j <= cnt && a[j] == a[i]) j++;
			now = 1ll * now * inv[j - i] % P;
		}
	}
	ans = (ans + 1ll * now * getdp()) % P;
}
void work(int lft, int lst, int depth) {
	if (lft == 0) {
		check(depth - 1);
		return;
	}
	if (lst == 0) return;
	work(lft, lst - 1, depth);
	if (lft >= lst) {
		a[depth++] = lst;
		work(lft - lst, lst, depth);
	}
}
int main() {
	read(n), read(m);
	fac[0] = inv[0] = 1;
	if (n > m) swap(n, m);
	for (int i = 1; i <= m; i++) {
		fac[i] = 1ll * fac[i - 1] * i % P;
		inv[i] = power(fac[i], P - 2);
	}
	two[0] = 1;
	for (int i = 1; i <= n * m; i++)
		two[i] = two[i - 1] * 2 % P;
	for (int i = 1; i <= m; i++)
	for (int j = 1; j <= m; j++)
		g[i][j] = __gcd(i, j);
	work(n, n, 1);
	writeln(1ll * ans * inv[n] % P * inv[m] % P);
	return 0;
}