結果

問題 No.310 2文字しりとり
ユーザー antaanta
提出日時 2015-12-09 14:45:53
言語 C++11
(gcc 11.4.0)
結果
AC  
実行時間 227 ms / 6,000 ms
コード長 8,529 bytes
コンパイル時間 1,113 ms
コンパイル使用メモリ 100,344 KB
実行使用メモリ 4,352 KB
最終ジャッジ日時 2023-10-13 08:40:43
合計ジャッジ時間 3,500 ms
ジャッジサーバーID
(参考情報)
judge12 / judge11
このコードへのチャレンジ(β)

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 1 ms
4,352 KB
testcase_01 AC 1 ms
4,348 KB
testcase_02 AC 1 ms
4,352 KB
testcase_03 AC 1 ms
4,352 KB
testcase_04 AC 2 ms
4,352 KB
testcase_05 AC 1 ms
4,348 KB
testcase_06 AC 2 ms
4,348 KB
testcase_07 AC 2 ms
4,352 KB
testcase_08 AC 2 ms
4,352 KB
testcase_09 AC 2 ms
4,348 KB
testcase_10 AC 1 ms
4,352 KB
testcase_11 AC 1 ms
4,348 KB
testcase_12 AC 2 ms
4,352 KB
testcase_13 AC 1 ms
4,348 KB
testcase_14 AC 2 ms
4,348 KB
testcase_15 AC 1 ms
4,348 KB
testcase_16 AC 1 ms
4,348 KB
testcase_17 AC 1 ms
4,348 KB
testcase_18 AC 2 ms
4,348 KB
testcase_19 AC 3 ms
4,348 KB
testcase_20 AC 4 ms
4,352 KB
testcase_21 AC 192 ms
4,348 KB
testcase_22 AC 227 ms
4,348 KB
testcase_23 AC 226 ms
4,348 KB
testcase_24 AC 225 ms
4,352 KB
testcase_25 AC 224 ms
4,348 KB
testcase_26 AC 46 ms
4,348 KB
testcase_27 AC 66 ms
4,352 KB
権限があれば一括ダウンロードができます

ソースコード

diff #

#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

typedef unsigned long long ull;

mint dot(const mint *a, const mint *b, int n) {
	const int K = 16;
	static_assert((ull)mint::Mod * mint::Mod < ~0ULL / (K + 1), "K is too large");
	ull sum = 0;
	int i;
	for(i = 0; i + K <= n; ) {
		rep(j, K) {
			sum += (ull)a[i].x * b[i].x;
			++ i;
		}
		sum %= mint::Mod;
	}
	for(; i < n; ++ i)
		sum += (ull)a[i].x * b[i].x;
	return mint((int)(sum % mint::Mod));
}

//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(vector<mint> s, vector<mint> &C) {
	int N = (int)s.size();
	reverse(s.begin(), s.end());	//逆順にしておく
	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 - 1 - n];
		if(L > 0) d += dot(&C[1], &s[N - 1 - n + 1], L);
		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].x = (C[m + i].x + (ull)coeff.x * B[i].x) % mint::Mod;
			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();
}

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);
	vector<ull> tmp(n);
	randomModIntVector(u);
	randomModIntVector(b);

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

	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) {
		tmp.assign(n, 0);

		ull sumb_t = 0;
		rep(i, n) {
			b[i] *= D[i];
			sumb_t += b[i].x;
		}

		if(src != -1) {
			if(dst < n && src < n)
				tmp[dst] += mint::Mod - b[src].x;
		}

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

		unsigned negsumb = mint::Mod - sumb_t % mint::Mod;
		rep(i, n)
			b[i].x = (tmp[i] + (ull)diag1[i].x * b[i].x + negsumb) % mint::Mod;
		
		uTAib[k] = dot(&u[0], &b[0], n);
	}

	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<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;

			edges[i] = {A, B};
			-- outDeg[A];
			-- inDeg[B];
		}

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

		//まず、eulerianかどうか判定する
		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(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;
}
0