C++11(clang++ 3.9) 解法, 执行用时: 284ms, 内存消耗: 1116K, 提交时间: 2019-07-06 08:31:21

#include<bits/stdc++.h> 
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rep1(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
#define rep2(i,l,u) for(int (i)=(int)(l);(i)<(int)(u);++(i))
using namespace std;

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;

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

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


void computeMinimumPolynomialForLinearlyRecurrentSequence(const vector<mint> &a, vector<mint> &phi) {
	int n2 = (int)a.size(), n = n2 / 2;
	int L = berlekampMessey(a, phi);
	reverse(phi.begin(), phi.begin() + (L + 1));
}

struct RandomModInt {
	default_random_engine re;
	uniform_int_distribution<int> dist;
	RandomModInt() : re(random_device {}()), dist(1, mint::Mod - 1) { }
	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);
	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) {
	vector<int> loops(N, 1);
	vector<int> outDeg(N, N), inDeg(N, N);
	if(src != -1) {
		++ outDeg[dst];
		++ inDeg[src];
	}
	for(auto e : edges) {
		-- outDeg[e.first];
		-- inDeg[e.second];
		if(e.first == e.second)
			-- loops[e.first];
	}

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

	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];
		res *= fact[deg - 1];
	}

	return res;
}

int main() {
	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");
		return 0;
	}

	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");
		return 0;
	}

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

	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);
	return 0;
}