伝説の良難問。
現在でこそ AC 数 30 人で解説記事も豊富にあるが、当時は AC 数 3 人という状況で解説も無い中で、必死に 1 週間かけて通した想い出の問題。
問題概要
正の整数 と 以上の整数 が与えられる。
{} から {} への写像 の組であって、以下の条件を満たすものが何通りあるか、1000000007 で割ったあまりを求めよ。
- id を恒等写像として、 = id
制約
考えたこと
順列に関する問題!!!
順列を写像として見て扱うのは、以下の問題でもお馴染み。
まず、 が全単射であることから、 は少なくとも全単射でなければならず、 については
の場合がコーナーケースとなっている。なお、コーナーケースの場合は の個数がイレギュラーであることを除けば、以下の議論はそう大差ない。
式変形
がともに全単射である場合 (順列となっている場合) について、式変形してみる。全単射なので逆写像が存在する。よって や を考えてよい。左から を、右から をかけることで、以下のようになる。
= id
⇔
ここで改めて とおくと、 の組と、 の組とは一対一に対応するので、結局
を満たすような順列 の組を求めればよいことになる。なお、 が負である場合には、改めて として、
となるので、 は 0 以上であるとしてよい (絶対値をとった値で考えてよい)。以上をまとめると、 として、記号を置き直して、
を満たすような順列 の組を求める問題ということになる。
数え上げの方針
ここまでくると簡単に思える方もいるらしい...僕はここまではむしろ一瞬だったのに、ここから詰め切るのに 1 週間かかった ( というコーナーケースには気づいてなかった)。
方針としては、各順列 に対して
を満たすような と をそれぞれ数え上げて、その積をとる。それを 通りの に対して総和をとる、という方針をとることにする。 の方が一般的なので について考える。
順列 f, h を巡回群の直積に分解して考える
整数問題では、整数を素因数分解できることから、まず素数の場合を解いて、素数べきの場合を解いて、最後に全体をまとめるようなアプローチをよくする。順列に対しても、巡回群の直積に分解できることから、まず巡回群の場合を解くのは自然だと思われる。
一旦仮に、 を位数 の巡回群であるとしてみよう ( の定義域のサイズは一時的に ではなく であるとする)。
このとき、たとえば と とが互いに素であったならば、 と とは一対一対応することがわかる (より正確にいえば、 に対して を返す写像は、「位数 の巡回群の集合」から「位数 の巡回群の集合」への全単射である。一般の場合には、
- の位数 は、 として、 となる
- と置くと と は互いに素であり、 を位数 の巡回群として考えている限りは、 を満たすような の組の個数を数えても等しい
ということがわかる。よって代わりにこれを数える。 が位数 の巡回群であるとき、 は「位数 な巡回群の 個の直積」となる。さて、まず、そのような を一つ固定したときに、それに対応する が何個あるのかを考えよう。
- は位数 の巡回群であるが、 を計算すると、 の各元は 個先の元へと移ることになる (つまり位数は となる)
- を固定したとき、 の各元が 個進んでどの元になるのかは一意に決まるので、よって において各元は 個のメンバをもつ 個の「数珠」に分かれることになる
- そこから を構成することは、各数珠の順序を固定した上で、数珠を繋ぎ直すことに相当する
- 個の数珠の順序を決める: 通り
- 個の数珠それぞれについて、中身の offset を決める: 通り
ということで、「位数 の巡回群の 個の直積」であるような に対して、「位数 の巡回群」であるような の個数は、 : に対応することがわかる。
DP へ
以上の考察から、 として、 を満たす の組の個数を求めるために、以下の DP を回すことを考えよう:
- dp[ ][ ] := を の順列としたときに、 の位数が 以下であるようなものの個数
そして、次のような更新式が立つ。 の部分は挿入 DP の考え方を使っていて、 の部分は「 個のものを グループに分けるグループ分けの個数」である。
- dp[ ] [ ] += dp[ ][ ] × × × × ( となるような の個数) × ( となるような の個数)
( を、位数 の巡回群の 個直積の一つとする)
最後の DP
最後に、 となるような の個数を考える。各 に対して、 は位数 の巡回群の 個の直積である。さて、 を構成する巡回群の位数を とするとき、
- として、 となる
- よって、 の倍数 のうち、 となるものを考えれば良い
そのような を とする。このとき、
とする。以下のような DP が立つ。各 に対して として、
と更新できる。 の部分は挿入 DP の思想に従って に新たに 個の巡回群 (位数 ) を付け加えることを表す (ここを としてしまい、バグに苦しんだ)。 の部分は、 をなす 個の巡回群を 個のグループに分ける方法の個数を表す。
コード
必死のデバッグの跡
#include <iostream> #include <vector> #include <algorithm> #include <map> using namespace std; long long GCD(long long x, long long y) { if (y == 0) return x; else return GCD(y, x % y); } // modint: mod 計算を int を扱うように扱える構造体 template<int MOD> struct Fp { long long val; constexpr Fp(long long v = 0) noexcept : val(v % MOD) { if (val < 0) val += MOD; } constexpr int getmod() { return MOD; } constexpr Fp operator - () const noexcept { return val ? MOD - val : 0; } constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; } constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; } constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; } constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; } constexpr Fp& operator += (const Fp& r) noexcept { val += r.val; if (val >= MOD) val -= MOD; return *this; } constexpr Fp& operator -= (const Fp& r) noexcept { val -= r.val; if (val < 0) val += MOD; return *this; } constexpr Fp& operator *= (const Fp& r) noexcept { val = val * r.val % MOD; return *this; } constexpr Fp& operator /= (const Fp& r) noexcept { long long a = r.val, b = MOD, u = 1, v = 0; while (b) { long long t = a / b; a -= t * b; swap(a, b); u -= t * v; swap(u, v); } val = val * u % MOD; if (val < 0) val += MOD; return *this; } constexpr bool operator == (const Fp& r) const noexcept { return this->val == r.val; } constexpr bool operator != (const Fp& r) const noexcept { return this->val != r.val; } friend constexpr ostream& operator << (ostream &os, const Fp<MOD>& x) noexcept { return os << x.val; } friend constexpr Fp<MOD> modpow(const Fp<MOD> &a, long long n) noexcept { if (n == 0) return 1; auto t = modpow(a, n / 2); t = t * t; if (n & 1) t = t * a; return t; } }; // 二項係数ライブラリ template<class T> struct BiCoef { vector<T> fact_, inv_, finv_; constexpr BiCoef() {} constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) { init(n); } constexpr void init(int n) noexcept { fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1); int MOD = fact_[0].getmod(); for(int i = 2; i < n; i++){ fact_[i] = fact_[i-1] * i; inv_[i] = -inv_[MOD%i] * (MOD/i); finv_[i] = finv_[i-1] * inv_[i]; } } constexpr T com(int n, int k) const noexcept { if (n < k || n < 0 || k < 0) return 0; return fact_[n] * finv_[k] * finv_[n-k]; } constexpr T fact(int n) const noexcept { if (n < 0) return 0; return fact_[n]; } constexpr T inv(int n) const noexcept { if (n < 0) return 0; return inv_[n]; } constexpr T finv(int n) const noexcept { if (n < 0) return 0; return finv_[n]; } }; const int MOD = 1000000007; using mint = Fp<MOD>; BiCoef<mint> bc; // e*k 人を、e 人の k グループに分ける方法の数 mint grouping(int e, int k) { return bc.fact(k*e) * bc.finv(k) * modpow(bc.finv(e), k); } // res[k] = h を位数 e の巡回群の k 個の直積からなる 1 つの順列としたときに // f^p = h となる f が何個あるか vector<mint> sub(int N, long long p, long long e) { int K = N / e; // f を構成する巡回群の位数 d として、e = d / GCD(d, p) となる d のみを考える vector<long long> D; for (long long d = e; d <= N; d += e) { if (e * GCD(d, p) == d) D.push_back(d); } // DP vector<vector<mint>> dp(D.size()+1, vector<mint>(K+1, 0)); dp[0][0] = 1; for (int i = 0; i < D.size(); ++i) { int d = D[i]; int g = d / e; mint mul = bc.fact(g - 1) * modpow(mint(e), g - 1); for (int k = 0; k <= K; ++k) { mint fac = 1; for (int l = 0; k + l * g <= K; ++l) { dp[i + 1][k + l * g] += dp[i][k] * bc.com(k + l * g, k) * grouping(g, l) * fac; fac *= mul; } } } return dp[D.size()]; } mint solve(int N, long long X, long long Y, long long Z) { if (X + Y + Z == 0) { auto gdp = sub(N, 2, 1); return gdp[N] * modpow(mint(N), N); } long long t = X - Y + Z; if (t < 0) t = -t; vector<vector<mint>> dp(N+1, vector<mint>(N+1, 0)); dp[0][0] = 1; for (int e = 1; e <= N; ++e) { const vector<mint>& fdp = sub(N, t, e); const vector<mint>& gdp = sub(N, 2, e); for (int n = 0; n <= N; ++n) { mint fac = 1; for (int k = 0; n + k * e <= N; ++k) { dp[e][n + k * e] += dp[e - 1][n] * bc.com(n + k * e, n) * grouping(e, k) * fac * fdp[k] * gdp[k]; fac *= bc.fact(e - 1); } } } return dp[N][N]; } int main() { bc.init(5100); int N; long long X, Y, Z; while (cin >> N >> X >> Y >> Z) cout << solve(N, X, Y, Z) << endl; }