#62272 (C++11) No.310 2文字しりとり

• 問題: No.310 2文字しりとり
• ユーザー: anta
• 提出日時: 2015-12-03 19:44:38
• 言語: C++11 (gcc 11.4.0)
• 結果: AC
• 実行時間: 439 ms / 6,000 ms
• コード長: 9,243 bytes
• コンパイル時間: 1,671 ms
• コンパイル使用メモリ: 106,180 KB
• 実行使用メモリ: 6,272 KB
• 最終ジャッジ日時: 2024-09-14 08:30:49
• 合計ジャッジ時間: 5,176 ms
• ジャッジサーバーID （参考情報）: judge2 / judge3

入力	結果	実行時間; 実行使用メモリ
testcase_00	AC	1 ms; 5,248 KB
testcase_01	AC	1 ms; 5,376 KB
testcase_02	AC	2 ms; 5,376 KB
testcase_03	AC	2 ms; 5,376 KB
testcase_04	AC	2 ms; 5,376 KB
testcase_05	AC	2 ms; 5,376 KB
testcase_06	AC	2 ms; 5,376 KB
testcase_07	AC	2 ms; 5,376 KB
testcase_08	AC	2 ms; 5,376 KB
testcase_09	AC	2 ms; 5,376 KB
testcase_10	AC	2 ms; 5,376 KB
testcase_11	AC	2 ms; 5,376 KB
testcase_12	AC	1 ms; 5,376 KB
testcase_13	AC	2 ms; 5,376 KB
testcase_14	AC	2 ms; 5,376 KB
testcase_15	AC	2 ms; 5,376 KB
testcase_16	AC	2 ms; 5,376 KB
testcase_17	AC	2 ms; 5,376 KB
testcase_18	AC	2 ms; 5,376 KB
testcase_19	AC	4 ms; 5,376 KB
testcase_20	AC	5 ms; 5,376 KB
testcase_21	AC	394 ms; 6,144 KB
testcase_22	AC	437 ms; 6,272 KB
testcase_23	AC	438 ms; 6,144 KB
testcase_24	AC	439 ms; 6,272 KB
testcase_25	AC	434 ms; 6,272 KB
testcase_26	AC	84 ms; 5,376 KB
testcase_27	AC	122 ms; 5,376 KB

コンパイルメッセージ

main.cpp: In function ‘int main()’:
main.cpp:368:26: warning: format ‘%d’ expects argument of type ‘int’, but argument 2 has type ‘mint’ {aka ‘ModInt<1000000007>’} [-Wformat=]
  368 |                 printf("%d\n", ans);
      |                         ~^     ~~~
      |                          |     |
      |                          int   mint {aka ModInt<1000000007>}
main.cpp:282:22: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  282 |                 scanf("%d%d", &N, &M);
      |                 ~~~~~^~~~~~~~~~~~~~~~
main.cpp:289:30: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  289 |                         scanf("%d%d", &A, &B), -- A, -- B;
      |                         ~~~~~^~~~~~~~~~~~~~~~

ソースコード

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <numeric>
#include <vector>
#include <random>
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <cassert>
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
#define reu(i,l,u) for(int (i)=(int)(l);(i)<(int)(u);++(i))
using namespace std;

////////////////////////////////////////////////////////
//ModInt

template<int MOD>
struct ModInt {
	static const int Mod = MOD;
	unsigned x;
	ModInt() : x(0) {}
	ModInt(signed sig) { int sigt = sig % MOD; if(sigt < 0) sigt += MOD; x = sigt; }
	ModInt(signed long long sig) { int sigt = sig % MOD; if(sigt < 0) sigt += MOD; x = sigt; }
	int get() const { return (int)x; }

	ModInt &operator+=(ModInt that) { if((x += that.x) >= MOD) x -= MOD; return *this; }
	ModInt &operator-=(ModInt that) { if((x += MOD - that.x) >= MOD) x -= MOD; return *this; }
	ModInt &operator*=(ModInt that) { x = (unsigned long long)x * that.x % MOD; return *this; }
	ModInt &operator/=(ModInt that) { return *this *= that.inverse(); }

	ModInt operator+(ModInt that) const { return ModInt(*this) += that; }
	ModInt operator-(ModInt that) const { return ModInt(*this) -= that; }
	ModInt operator*(ModInt that) const { return ModInt(*this) *= that; }
	ModInt operator/(ModInt that) const { return ModInt(*this) /= that; }

	ModInt inverse() const {
		signed a = x, b = MOD, u = 1, v = 0;
		while(b) {
			signed t = a / b;
			a -= t * b; std::swap(a, b);
			u -= t * v; std::swap(u, v);
		}
		if(u < 0) u += Mod;
		ModInt res; res.x = (unsigned)u;
		return res;
	}

	bool operator==(ModInt that) const { return x == that.x; }
	bool operator!=(ModInt that) const { return x != that.x; }
	ModInt operator-() const { ModInt t; t.x = x == 0 ? 0 : Mod - x; return t; }
};
typedef ModInt<1000000007> mint;

////////////////////////////////////////////////////////
//Black box linear algebra

//Berlekamp-Massey algorithm
//(\sum_{j=0}^L C[j] s[i + L - j]) = 0 for all 0 <= i < N - L
//となるような最小の L, C を返す。C[0] = 1 である
int berlekampMessey(const vector<mint> &s, vector<mint> &C) {
	int N = (int)s.size();
	C.assign(N + 1, mint());
	vector<mint> B(N + 1, mint());
	C[0] = B[0] = 1;
	int degB = 0;
	vector<mint> T;
	int L = 0, m = 1;
	mint b = 1;
	for(int n = 0; n < N; ++ n) {
		mint d = s[n];
		for(int i = 1; i <= L; ++ i)
			d += C[i] * s[n - i];
		if(d == mint()) {
			++ m;
		} else {
			if(2 * L <= n)
				T.assign(C.begin(), C.begin() + (L + 1));
			mint coeff = -d * b.inverse();
			for(int i = 0; i <= degB; ++ i)
				C[m + i] += coeff * B[i];
			if(2 * L <= n) {
				L = n + 1 - L;
				B.swap(T);
				degB = (int)B.size() - 1;
				b = d;
				m = 1;
			} else {
				++ m;
			}
		}
	}
	C.resize(L + 1);
	return L;
}

//体上の数列のminimum polynomialを計算する。
//(\sum_{j=0}^d phi[j] s[i + j]) = 0 for all i
void computeMinimumPolynomialForLinearlyRecurrentSequence(const vector<mint> &a, vector<mint> &phi) {
	int n2 = (int)a.size(), n = n2 / 2;
	assert(n2 % 2 == 0);
	int L = berlekampMessey(a, phi);
	reverse(phi.begin(), phi.begin() + (L + 1));
}

struct RandomModInt {
	default_random_engine re;
	uniform_int_distribution<int> dist;
#ifndef _DEBUG
	RandomModInt() : re(random_device{}()), dist(1, mint::Mod - 1) { }
#else
	RandomModInt() : re(), dist(1, mint::Mod - 1) { }
#endif
	mint operator()() {
		mint r;
		r.x = dist(re);
		return r;
	}
} randomModInt;

void randomModIntVector(vector<mint> &v) {
	int n = (int)v.size();
	for(int i = 0; i < n; ++ i)
		v[i] = randomModInt();
}


////////////////////////////////////////////////////////
//ユーティリティー
struct UnionFind {
	vector<int> data;
	void init(int n) { data.assign(n, -1); }
	bool unionSet(int x, int y) {
		x = root(x); y = root(y);
		if(x != y) {
			if(data[y] < data[x]) swap(x, y);
			data[x] += data[y]; data[y] = x;
		}
		return x != y;
	}
	bool findSet(int x, int y) { return root(x) == root(y); }
	int root(int x) { return data[x] < 0 ? x : data[x] = root(data[x]); }
	int size(int x) { return -data[root(x)]; }
};

////////////////////////////////////////////////////////
//解答メイン部分

#ifdef MY_LOCAL_RUN
#define RANDOM_TEST
#define NAIVE_CHECK
//#define TIME_CHECK
#endif

mint dot(const vector<mint> &a, const vector<mint> &b) {
	int n = (int)a.size();
	assert(b.size() == n);
	mint sum;
	rep(i, n)
		sum += a[i] * b[i];
	return sum;
}

mint computeDeterminant(int N, const vector<mint> &diag, const vector<pair<int,int> > &validEdges, int src, int dst) {
	int n = N - 1;
	if(n == 0)
		return 1;
	vector<mint> D(n);
	vector<mint> m;
	randomModIntVector(D);
	vector<mint> u(n), b(n), res;
	randomModIntVector(u);
	randomModIntVector(b);

	vector<mint> uTAib(n * 2);
	uTAib[0] = dot(u, b);

	vector<mint> diag1(n);
	rep(i, n)
		diag1[i] = diag[i] + 1;

	vector<pair<short, short> > edges(validEdges.begin(), validEdges.end());

	for(int k = 1; k < n * 2; ++ k) {
		res.assign(n, mint());
		rep(i, n)
			b[i] *= D[i];

		mint sumb;
		rep(i, n)
			sumb += b[i];

		if(src != -1) {
			if(dst < n && src < n)
				res[dst] -= b[src];
		}

		for(auto e : edges)
			res[e.first] += b[e.second];

		rep(i, n)
			res[i] += diag1[i] * b[i] - sumb;

		b.swap(res);
		uTAib[k] = dot(u, b);
	}

	computeMinimumPolynomialForLinearlyRecurrentSequence(uTAib, m);
	//運が悪い場合は気にしないことにする
	//m = char(AD)
	assert(m.size() == n + 1);
	assert(m.back() == mint(1));
	mint detD = 1;
	for(int i = 0; i < n; ++ i)
		detD *= D[i];
	mint invdetD = detD.inverse();
	mint detA = m[0] * invdetD;
	if(n % 2 == 1)
		detA = mint() - detA;
	return detA;
}

mint countEulerianCyclesOnComplementGraph(int N, const vector<pair<int,int> > &edges, int src, int dst) {
	assert(N > 0);

	vector<int> loops(N, 1);
	vector<int> outDeg(N, N), inDeg(N, N);
	//サイクルにならない場合は辺を1つ足してサイクルにするようにする
	if(src != -1) {
		assert(src != dst);
		++ outDeg[dst];
		++ inDeg[src];
	}
	for(auto e : edges) {
		-- outDeg[e.first];
		-- inDeg[e.second];
		if(e.first == e.second)
			-- loops[e.first];
	}

	rep(i, N)
		assert(inDeg[i] == outDeg[i]);

	//BEST theoremによりオイラー閉路を数えられる
	//<https://en.wikipedia.org/wiki/BEST_theorem>

	vector<mint> fact(N + 1);
	fact[0] = 1;
	rer(n, 1, N)
		fact[n] = fact[n - 1] * n;

	//matrix tree theorem によって有向木を数える
	//<https://en.wikipedia.org/wiki/Kirchhoff%27s_theorem#Kirchhoff.27s_theorem_for_directed_multigraphs>
	//頂点 N - 1 を根とする。
	//1行削除されるため、それに接する辺も一緒に削除しておく
	//また、ループも削除しておく
	vector<pair<int, int> > validEdges;
	for(auto e : edges) if(e.first < N - 1 && e.second < N - 1 && e.first != e.second)
		validEdges.push_back(e);
	//なんとなくソートしておく
	sort(validEdges.begin(), validEdges.end());
	//行列の対角要素
	vector<mint> diag(N - 1);
	rep(i, N - 1)
		diag[i] = outDeg[i] - loops[i];

	mint res = computeDeterminant(N, diag, validEdges, src, dst);

	rep(i, N) {
		int deg = outDeg[i];
		assert(deg > 0);
		res *= fact[deg - 1];
	}

	return res;
}

int main() {
	do {
		int N, M;
		scanf("%d%d", &N, &M);
		vector<vector<bool> > g(N, vector<bool>(N, true));
		vector<int> outDeg(N, N), inDeg(N, N);
		vector<pair<int, int> > edges(M);

		rep(i, M) {
			int A, B;
			scanf("%d%d", &A, &B), -- A, -- B;

			assert(g[A][B]);
			g[A][B] = false;
			edges[i] = {A, B};
			-- outDeg[A];
			-- inDeg[B];
		}

		//空グラフは場合分けしておく
		if(M == N * N) {
			puts("1");
			continue;
		}

		//ちなみに、ある補グラフ上でのDFSやBFSは線形時間で可能である
		UnionFind uf; uf.init(N);
		rep(i, N) rep(j, N) if(g[i][j])
			uf.unionSet(i, j);

		//まず、eulerianかどうか判定する
		int mainCC = -1;
		int src = -1, dst = -1;
		bool eulerian = true;
		rep(i, N) {
			if(outDeg[i] == inDeg[i]) {
			} else if(outDeg[i] == inDeg[i] + 1) {
				eulerian &= src == -1;
				src = i;
			} else if(outDeg[i] == inDeg[i] - 1) {
				eulerian &= dst == -1;
				dst = i;
			} else {
				eulerian = false;
			}
			if(mainCC == -1 && outDeg[i] > 0)
				mainCC = i;
		}
		assert(mainCC != -1);

		vector<int> vertices;
		rep(i, N) {
			if(outDeg[i] > 0 || inDeg[i] > 0)
				eulerian &= uf.findSet(mainCC, i);
		}

		if(src != -1 || dst != -1)
			eulerian &= src != -1 && dst != -1;

		if(!eulerian) {
			puts("0");
			continue;
		}

		//孤立点を取り除いたグラフを作る
		vector<int> vertexIndex(N, -1);
		int V = 0;
		rep(i, N)
			if(outDeg[i] > 0 || inDeg[i] > 0)
				vertexIndex[i] = V ++;
		if(src != -1) {
			src = vertexIndex[src];
			dst = vertexIndex[dst];
			assert(src != -1 && dst != -1);
		}

		for(auto &e : edges) {
			e.first = vertexIndex[e.first];
			e.second = vertexIndex[e.second];
			if(e.first == -1 || e.second == -1)
				e = {-1, -1};
		}
		edges.erase(remove(edges.begin(), edges.end(), make_pair(-1, -1)), edges.end());

		mint ans = countEulerianCyclesOnComplementGraph(V, edges, src, dst);
		//サイクルの場合、どの辺から始めるかの分をかける
		if(src == -1)
			ans *= N * N - M;

		printf("%d\n", ans);

	} while(0);
	return 0;
}