////////////////////////////////////////////////////////
//テストとかジェネレーターとか全部残してある
#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

void truncateLeadingZerosOfPolynomial(const vector<mint> &p, vector<mint> &res) {
	int n = (int)p.size();
	while(n > 0 && p[n-1].x == 0) -- n;
	res.assign(p.begin(), p.begin() + n);
}

int extendedEuclideanAlgorithm(const vector<mint> &f, const vector<mint> &g, int k, vector<mint> &out_r, vector<mint> &out_t, bool traditional = false) {
	if((!f.empty() && f.back().x == 0) || (!g.empty() && g.back().x == 0)) {
		vector<mint> nf, ng;
		truncateLeadingZerosOfPolynomial(f, nf);
		truncateLeadingZerosOfPolynomial(g, ng);
		return extendedEuclideanAlgorithm(nf, ng, k, out_r, out_t, traditional);
	}
	if(f.size() < g.size())
		return extendedEuclideanAlgorithm(g, f, k, out_r, out_t, traditional);
	assert(k >= 0);

	int n[2];
	vector<mint> r[2], t[2];//, s[2];
	int nf = (int)f.size() - 1, ng = (int)g.size() - 1;
	n[0] = nf, n[1] = ng;
	mint ip0 = traditional || f.empty() ? mint(1) : f.back().inverse();
	mint ip1 = traditional || g.empty() ? mint(1) : g.back().inverse();
	int bufsize = n[0] == -1 ? 1 : n[0] + 1;
	rep(i, 2) {
		r[i].resize(bufsize);
//		s[i].resize(bufsize);
		t[i].resize(bufsize);
	}
	rer(j, 0, n[0]) r[0][j] = f[j] * ip0;
	rer(j, 0, n[1]) r[1][j] = g[j] * ip1;
//	s[0][0] = ip0;
	t[1][0] = ip1;

	if(n[0] < k) {
		out_r.swap(r[0]); out_r.resize(n[0] + 1);
		out_t.swap(t[0]); out_t.resize(1);
		return 0;
	}else if(n[1] < k) {
		out_r.swap(r[1]); out_r.resize(n[1] + 1);
		out_t.swap(t[1]); out_t.resize(1);
		return 1;
	}

	for(int i = 1; n[1] >= 0; ++ i) {
		swap(n[0], n[1]);
		r[0].swap(r[1]);
//		s[0].swap(s[1]);
		t[0].swap(t[1]);
		//r[0] = r_i; r[1] = r_{i-1}, r_{i+1}

		mint ilc = r[0][n[0]].inverse();
		assert(traditional || ilc.x == 1);
		//q_i     <- r_{i-1} quo r_i
		//r_{i+1} <- r_{i-1} - q_i r_i
		//s_{i+1} <- s_{i-1} - q_i s_i
		//t_{i+1} <- t_{i-1} - q_i t_i
		int nr = n[0], ns = ng - n[1], nt = nf - n[1];

		for(int j = n[1] - nr; j >= 0; -- j) {
			mint c = r[1][nr + j] * ilc;
			for(int l = 0; l <= nr; ++ l)
				r[1][j + l] -= c * r[0][l];
//			for(int l = 0; l <= ns; ++ l)
//				s[1][j + l] -= c * s[0][l];
			for(int l = 0; l <= nt; ++ l)
				t[1][j + l] -= c * t[0][l];
		}
		while(n[1] >= 0 && r[1][n[1]].x == 0)
			-- n[1];
		mint ip = n[1] == -1 ? mint(1) : r[1][n[1]].inverse();
		if(!traditional) {
			rer(j, 0, n[1]) r[1][j] *= ip;
//			rer(j, 0, ng - n[0]) s[1][j] *= ip;
			rer(j, 0, nf - n[0]) t[1][j] *= ip;
		}

		if(n[1] < k) {
			out_r.swap(r[1]); out_r.resize(n[1] + 1);
			out_t.swap(t[1]); out_t.resize(nf - n[0] + 1);
			return i+1;
		}
	}
	assert(false);
	return -1;
}

//体上の数列のminimum polynomialを計算する。
//入力がlinearly recurrentでない場合は、res.back() == 0 になる。

void computeMinimumPolynomialForLinearlyRecurrentSequence(const vector<mint> &a, vector<mint> &res) {
	int n2 = (int)a.size(), n = n2 / 2;
	assert(n2 % 2 == 0);
	vector<mint> x(n2 + 1), s, t;
	x[n2] = 1;
	int i = extendedEuclideanAlgorithm(x, a, n, s, t);
	if(s.size() == n2 + 1) {	//a is the zero sequence
		res.assign(1, mint());
		res[0] = mint(1);
		return;
	}
	
	int d = max((int)s.size(), (int)t.size() - 1);
	int e = 0;
	while(t[e].x == 0) ++ e;
	mint c = t[e].inverse();
	res.assign(d + 1, mint());
	reu(i, e, t.size())
		res[d - i] = t[i] * c;
}

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

bool verifyMinimalPolynomialForKrylovSubspace(const vector<vector<mint> > &Aib, const vector<mint> &m) {
	int n = (int)Aib[0].size();
	assert(!m.empty());
	int d = (int)m.size() - 1;
	assert(d < (int)Aib.size());
	vector<mint> sum(n);
	for(int i = 0; i <= d; ++ i) {
		mint coeff = m[i];
		const vector<mint> &v = Aib[i];
		for(int j = 0; j < n; ++ j)
			sum[j] += coeff * v[j];
	}
	for(int j = 0; j < n; ++ j)
		if(sum[j] != mint())
			return false;
	return true;
}

//誤った答えを返すことがある (ただし minimal polynomial の約数である)
//確実にしたい場合、verifyMinimalPolynomialForKrylovSubspaceでverifyされるまでイテレーションすればよい。
void minimalPolynomialForKrylovSubspace(const vector<vector<mint> > &Aib, vector<mint> &res) {
	int n = (int)Aib[0].size();
	assert(Aib.size() == n * 2);
	vector<mint> u(n);
	randomModIntVector(u);
	vector<mint> uTAib(n * 2);
	for(int i = 0; i < n * 2; ++ i)
		uTAib[i] = inner_product(u.begin(), u.end(), Aib[i].begin(), mint());
	computeMinimumPolynomialForLinearlyRecurrentSequence(uTAib, res);
}

//誤った答えを返すことがある (ただし minimal polynomial の約数である)
//A(b,res): res += Ab
template<typename BlackBox>
void blackBoxMinimalPolynomial(const BlackBox &A, int n, vector<mint> &m) {
	if(n == 0) {
		m.assign(1, mint(1));
		return;
	}
	vector<vector<mint> > Aib(n * 2, vector<mint>(n));
	vector<mint> &b = Aib[0];
	randomModIntVector(b);
	for(int i = 0; i < n * 2 - 1; ++ i)
		A(Aib[i], Aib[i + 1]);
	minimalPolynomialForKrylovSubspace(Aib, m);
}

template<typename BlackBox>
struct DiagonalPreconditionBox {
	const BlackBox &A;
	const vector<mint> D;

	DiagonalPreconditionBox(const BlackBox &A, const vector<mint> &D): A(A), D(D) { }

	//res = ADb
	void operator()(const vector<mint> &b, vector<mint> &res) const {
		int n = (int)b.size();
		assert(D.size() == n);
		vector<mint> Db(n);
		for(int i = 0; i < n; ++ i)
			Db[i] = D[i] * b[i];
		A(Db, res);
	}
};

//常に正しい答えを返す
template<typename BlackBox>
mint blackBoxDeterminant(const BlackBox &A, int n) {
	vector<mint> D(n);
	vector<mint> m;
	do {
		randomModIntVector(D);
		DiagonalPreconditionBox<BlackBox> AD(A, D);
		blackBoxMinimalPolynomial(AD, n, m);
		assert(!m.empty());
		if(m[0] == mint()) return mint();
	} while((int)m.size() < n + 1);
	//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;
}

////////////////////////////////////////////////////////
//ユーティリティー

#ifdef NDEBUG
#error assert is disabled!
#endif

template<typename T>
T readNatural(T lo, T up) {
	assert(0 <= lo && lo <= up);
	T x = 0;
	while(1) {
		int d = getchar();
		if(!('0' <= d && d <= '9')) {
			ungetc(d, stdin);
			break;
		}
		d -= '0';
		assert(d <= up && x <= (up - d) / 10);
		x = x * 10 + d;
	}
	assert(lo <= x && x <= up);
	return x;
}

void readSpace() { int c = getchar(); assert(c == ' '); }
static bool read_eof = false;
void readEOL() {
	int c = getchar();
	if(c == EOF) {
		assert(!read_eof);
		read_eof = true;
	}else {
		assert(c == '\n');
	}
}
void readEOF() {
	assert(!read_eof);
	int c = getchar();
	assert(c == EOF);
	read_eof = true;
}

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 TIME_CHECK
//#define NAIVE_CHECK
#endif

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 = blackBoxDeterminant([N, &diag, &validEdges, src, dst](const vector<mint> &b, vector<mint> &res) {
		assert(b.size() == N - 1 && res.size() == N - 1);
		//A[i][j] where i != j:
		//= -(1 + (i == dst && j == src) - ((i,j) in edges))
		mint sumb;
		rep(i, N - 1)
			sumb += b[i];
		//自分以外の場所から1引かれる
		rep(i, N - 1)
			res[i] -= sumb - b[i];
		if(src != -1) {
			//dst -> src の辺
			if(dst < N - 1 && src < N - 1)
				res[dst] -= b[src];
		}
		for(auto e : validEdges)
			res[e.first] += b[e.second];

		//A[i][i] = diag[i]
		rep(i, N - 1)
			res[i] += diag[i] * b[i];
	}, N - 1);

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

	return res;
}


mint naiveSolve(const vector<vector<bool> > &g) {
	int N = (int)g.size();
	vector<pair<int, int> > aliveEdges;
	rep(i, N) rep(j, N) if(g[i][j])
		aliveEdges.push_back({i, j});
	int E = (int)aliveEdges.size();
	assert(E <= 25);
	if(E == 0)
		return 1;
	vector<mint> dp(N << E);
	rep(i, N)
		dp[0 * N + i] += 1;
	rep(s, 1 << E) rep(i, N) {
		mint x = dp[s * N + i];
		if(x.x == 0) continue;
		rep(j, E) if((~s >> j & 1) && aliveEdges[j].first == i)
			dp[(s | 1 << j) * N + aliveEdges[j].second] += x;
	}
	mint res;
	rep(i, N)
		res += dp[((1 << E)-1) * N + i];
	return res;
}


void naiveCheck(const vector<vector<bool> > &g, mint ans) {
#ifdef NAIVE_CHECK
	mint ansnaive = naiveSolve(g);
	if(ansnaive.get() != ans.get()) {
		cerr << ansnaive.get() << " != " << ans.get() << endl;
		abort();
	}
#endif
}

#ifdef TIME_CHECK
struct FILETIME {
	unsigned dwLowDateTime, dwHighDateTime;
};
extern "C" void* __stdcall GetCurrentProcess(void);
extern "C" int __stdcall GetProcessTimes(void *hProcess, FILETIME *lpCreationTime, FILETIME *lpExitTime, FILETIME *lpKernelTime, FILETIME *lpUserTime);
void getCPUTime(double &userTime, double &systemTime) {
	void *handle = GetCurrentProcess();
	FILETIME dummy1, dummy2, kernel, user;
	GetProcessTimes(handle, &dummy1, &dummy2, &kernel, &user);
	userTime = user.dwHighDateTime * 429.4967296 + user.dwLowDateTime * 1e-7;
	systemTime = kernel.dwHighDateTime * 429.4967296 + kernel.dwLowDateTime * 1e-7;
}
double getCPUTime() {
	double user, sys;
	getCPUTime(user, sys);
	return user + sys;
}
struct CPUTimeIt {
	double user, sys;
	const char *msg;

	CPUTimeIt(const char *msg_): msg(msg_) { getCPUTime(user, sys); }
	~CPUTimeIt() {
		double userEnd, sysEnd;
		getCPUTime(userEnd, sysEnd);
		fprintf(stderr, "%s: user %.6fs / sys %.6fs\n", msg, userEnd - user, sysEnd - sys);
	}

	operator bool() { return false; }
};
#define CPUTIMEIT(s) if(CPUTimeIt cputimeit_##__LINE__ = s); else
#endif

namespace test_det {
void test();
void test_classical_time();
}

int main() {
//	freopen("C:/test/sagyo/874/test_0.txt", "r", stdin);
#ifdef RANDOM_TEST
	default_random_engine re(random_device{}());
#endif

#ifdef TIME_CHECK
	test_det::test_classical_time();
#endif

	do
#ifdef TIME_CHECK
	CPUTIMEIT("main")
#endif
	{
		const int MaxN = 4000, MaxM = 4000;
#ifndef RANDOM_TEST
		int N = readNatural(1, MaxN); readSpace();
		int M = readNatural(0, MaxM); readEOL();
#else
#ifndef TIME_CHECK
		int N = uniform_int_distribution<int>(3, 3)(re);
		int M = uniform_int_distribution<int>(0, N * N)(re);
#else
		int N = MaxN;
		int M = MaxM;
#endif
#endif
		vector<vector<bool> > g(N, vector<bool>(N, true));
		vector<int> outDeg(N, N), inDeg(N, N);
		vector<pair<int, int> > edges(M);
#if defined(RANDOM_TEST) && defined(TIME_CHECK)
		vector<int> randPerm(N);
		iota(randPerm.begin(), randPerm.end(), 0);
		shuffle(randPerm.begin(), randPerm.end(), re);
#endif
		rep(i, M) {
#ifndef RANDOM_TEST
			int A = readNatural(1, N) - 1; readSpace();
			int B = readNatural(1, N) - 1; readEOL();
#else
#ifndef TIME_CHECK
			int A, B;
			do {
				A = uniform_int_distribution<int>(0, N - 1)(re);
				B = uniform_int_distribution<int>(0, N - 1)(re);
			} while(!g[A][B]);
#else
			int A = randPerm[i % N], B = randPerm[(i + 1) % N];
#endif
#endif

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

		//空グラフは場合分けしておく
		if(M == N * N) {
			naiveCheck(g, 1);
			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) {
			naiveCheck(g, 0);
			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;

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

#ifdef MY_LOCAL_RUN
	} while(1);
#else
	} while(0);
#endif

	readEOF();
	return 0;
}

////////////////////////////////////////////////////////
//テスト

namespace test_det {

namespace classical {

typedef mint Num;
typedef vector<Num> Vec;
typedef vector<Vec> Mat;
 
Num determinant(Mat A) {
	int n = (int)A.size();
	if(!A.empty() && n != A[0].size()) return Num();
	Num det = 1;
	for(int i = 0; i < n; ++ i) {
		int pi = i;
		while(pi < n && A[pi][i] == Num()) ++ pi;
		if(pi == n) return Num();
		A[pi].swap(A[i]);
		det *= A[i][i];
		if(pi != i) det = -det;
		Num inv = Num(1) / A[i][i];
		for(int j = i + 1; j < n; ++ j) {
			Num t = A[j][i] * inv;
			if(t == Num()) continue;
			for(int k = i; k < n; ++ k)
				A[j][k] -= t * A[i][k];
		}
	}
	return det;
}

}

void test() {
	default_random_engine re;
	uniform_int_distribution<int> dist(1, mint::Mod);
	rep(ii, 100000) {
		const int X = mint::Mod - 1;
		int n = uniform_int_distribution<int>(1, 10)(re);
		classical::Mat A(n, classical::Vec(n));
		rep(i, n) rep(j, n)
			A[i][j].x = uniform_int_distribution<int>(0, X)(re);
		mint det1 = classical::determinant(A);
		mint det2 = blackBoxDeterminant([n, &A](const vector<mint> &b, vector<mint> &res) {
			assert(res.size() == n && count(res.begin(), res.end(), mint()) == n);
			rep(i, n)
				res[i] = inner_product(b.begin(), b.end(), A[i].begin(), res[i]);
		}, n);
		if(det1 != det2) {
			cerr << det1.get() << " != " << det2.get() << endl;
			rep(i, n) {
				rep(j, n)
					cerr << A[i][j].get() << " ";
				cerr << endl;
			}
		}
	}
}

void test_classical_time() {
#ifdef TIME_CHECK
	default_random_engine re{777};
	int N = 1000;
	cerr << "N = " << N << endl;
	vector<vector<mint> > A(N, vector<mint>(N, mint(1)));
	rep(i, N)
		A[i][i] = mint();
	rep(k, N) {
		int i = uniform_int_distribution<int>(0, N - 1)(re);
		int j = uniform_int_distribution<int>(0, N - 1)(re);
		A[i][j] = randomModInt();
	}
	mint det;
	CPUTIMEIT("det") {
		det = classical::determinant(A);
	}
	cerr << "det = " << det.get() << endl;
#else
	abort();
#endif
}

}

////////////////////////////////////////////////////////
//テストケースジェネレーター
#if 0
//テストケース生成のためのもの

#include <string>
#include <vector>
#include <algorithm>
#include <numeric>
#include <set>
#include <map>
#include <iostream>
#include <sstream>
#include <fstream>
#include <limits>
#include <functional>
#include <random>
#include <tuple>
#include <memory>
#include <cstdio>
#include <cassert>
using namespace std;

#ifdef NDEBUG
#undef NDEBUG
#endif

#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))

#include "C:\Dropbox\backup\implements\Util\TimeIt_Simple.hpp"

//問題の識別子
const string problem_name = "874";
//テストケースを生成するディレクトリ
const string testcase_dir = "C:\\test\\sagyo\\" + problem_name;
//テストケースの拡張子。この拡張子のファイルはなんでも自動的に消されるので注意
const string testcase_postfix = ".txt";

struct RNG {
	default_random_engine re;
	RNG(): re() {}

	void addSeed(unsigned x) {
		re.seed(re() ^ x);
	}

	unsigned operator()() {
		return re();
	}
	int operator()(int N) {
		return uniform_int_distribution<int>(0, N - 1)(re);
	}
	int operator()(int min, int max) {
		return uniform_int_distribution<int>(min, max)(re);
	}
};

class TestCase {
public:
	TestCase(int no = -1) : no(no) {}

	typedef vector<pair<int, int> > Edges;
	static const int MaxN = 4000, MaxM = 4000;

	virtual void generate(RNG &rng) = 0;

	virtual string getName() const {
		return "test_" + to_string(no);
	}

	void write(ostream &os) const {
		int N;
		vector<pair<int, int> > edges;
		tie(N, edges) = get();
		int M = (int)edges.size();
		assert(1 <= N && N <= MaxN);
		assert(0 <= M && M <= MaxM);
		os << N << ' ' << M << '\n';
		vector<vector<bool> > g(N, vector<bool>(N, true));
		for(auto p : edges) {
			assert(0 <= p.first && p.first < N);
			assert(0 <= p.second && p.second < N);
			assert(g[p.first][p.second]);
			g[p.first][p.second] = false;
			os << p.first + 1 << ' ' << p.second + 1 << '\n';
		}
		os.flush();
	}

protected:
	virtual tuple<int,Edges> get() const = 0;

private:
	int no;
};

class Manual : public TestCase {
protected:
	int N;
	Edges edges;

public:
	Manual(int no, int N, Edges edges) : TestCase(no), N(N), edges(edges) {}

	void generate(RNG &rng) override {
		//nope
	}

protected:
	tuple<int,Edges> get() const override { return make_tuple(N, edges); }
};

class RandomBalanced : public TestCase {
protected:
	int N, M;
	bool notCycle;
	Edges edges;

public:
	RandomBalanced(int no, int N, int M, bool notCycle = false)
		: TestCase(no), N(N), M(M), notCycle(notCycle) {}

	void generate(RNG &rng) override {
		edges.clear();
		int E = 0;
		vector<int> v;
		vector<vector<bool> > g(N, vector<bool>(N, true));
		if(notCycle) ++ M;
		while(E < M) {
			int len, s;
			while(1) {
				len = rng(1, M - E);
				s = rng(N);
				v.clear();
				if(len == 1) {
					if(!g[s][s]) {
						continue;
					} else {
						v.push_back(s);
						break;
					}
				}
				v.push_back(s);
				bool ok = gen_dfs(1, rng, g, len, v);
				if(!ok)
					continue;
				break;
			}
			for(int k = 0; k < (int)v.size() - 1; ++ k)
				edges.push_back({v[k], v[k + 1]});
			g[v.back()][v.front()] = false;
			edges.push_back({v.back(), v.front()});
			E = (int)edges.size();
		}
		if(notCycle) {
			-- M;
			edges.erase(edges.begin() + rng((int)edges.size()));
		}
	}

private:
	bool gen_dfs(int k, RNG &rng, vector<vector<bool> > &g, int len, vector<int> &v) {
		if(k == len)
			return true;
		for(int trys = 0; trys < 100; ++ trys) {
			int j = rng(N);
			if(!g[v.back()][j])
				continue;
			g[v.back()][j] = false;
			if(k == len - 1 && !g[j][v.front()]) {
				g[v.back()][j] = true;
				continue;
			}
			v.push_back(j);
			if(gen_dfs(k + 1, rng, g, len, v))
				return true;
			v.pop_back();
			g[v.back()][j] = true;
		}
		return false;
	}

protected:
	tuple<int,Edges> get() const override { return make_tuple(N, edges); }
};

class RandomBigCycle : public TestCase {
protected:
	int N;
	Edges edges;

public:
	RandomBigCycle(int no, int N) : TestCase(no), N(N) {}

	void generate(RNG &rng) override {
		vector<int> perm(N);
		iota(perm.begin(), perm.end(), 0);
		shuffle(perm.begin(), perm.end(), rng.re);
		edges.clear();
		rep(i, N)
			edges.push_back({perm[i], perm[(i + 1) % N]});
	}

protected:
	tuple<int,Edges> get() const override { return make_tuple(N, edges); }
};


class RandomBidirected : public TestCase {
protected:
	int N;
	int halfM;
	Edges edges;

public:
	RandomBidirected(int no, int N, int halfM) : TestCase(no), N(N), halfM(halfM) {}

	void generate(RNG &rng) override {
		edges.clear();
		vector<vector<bool> > g(N, vector<bool>(N, true));
		for(int i = 0; i < halfM; ++ i) {
			int a, b;
			do {
				a = rng(N), b = rng(N);
			}while(a == b || !g[a][b]);
			edges.push_back({a, b});
			edges.push_back({b, a});
			g[a][b] = g[b][a] = false;
		}
	}

protected:
	tuple<int,Edges> get() const override { return make_tuple(N, edges); }
};


int main() {
	vector<unique_ptr<TestCase> > testcases;
	auto A = [&](TestCase *t) { testcases.emplace_back(t); };

	const int MaxN = TestCase::MaxN, MaxM = TestCase::MaxM;
	int no = 0;
	//サンプル
	A(new Manual(no ++, 2, {{0, 1}}));
	A(new Manual(no ++, 2, {}));
	A(new Manual(no ++, 2, {{1, 1}, {0, 1}, {0, 0}, {1, 0}}));
	A(new Manual(no ++, 3, {{0, 2}, {2, 0}, {2, 1}}));
	A(new Manual(no ++, 3, {}));
	A(new Manual(no ++, 10, {}));
	//コーナーケースたち
	A(new Manual(no ++, 1, {}));
	A(new Manual(no ++, 1, {{0, 0}}));
	A(new Manual(no ++, 2, {{0, 1}, {1, 0}}));	//未連結
	A(new Manual(no ++, 2, {{0, 0}, {0, 1}, {1, 0}}));	//孤立点がある
	A(new Manual(no ++, 4, {{0, 1}, {0, 2}}));	//unbalanced
	//ランダム
	A(new RandomBidirected(no ++, 5, 5));
	A(new RandomBalanced(no ++, 5, 10));
	A(new RandomBalanced(no ++, 5, 10, true));
	A(new RandomBalanced(no ++, 5, 15));
	A(new RandomBalanced(no ++, 5, 15, true));
	A(new RandomBidirected(no ++, 20, 50));
	A(new RandomBalanced(no ++, 20, 100));
	A(new RandomBalanced(no ++, 100, 200));
	A(new RandomBalanced(no ++, 70, MaxM));
	A(new RandomBalanced(no ++, 100, MaxM));
	//最大ケース1
	A(new Manual(no ++, MaxN, {}));
	//最大ケース2
	A(new RandomBidirected(no ++, MaxN, MaxM / 2));
	A(new RandomBalanced(no ++, MaxN, MaxM));
	A(new RandomBalanced(no ++, MaxN, MaxM, true));
	//最大ケース3
	A(new RandomBigCycle(no ++, MaxN));
	
	if(testcase_dir[0] != 'C') abort();
	if(testcase_postfix[0] != '.') abort();

	if(0) {
		cerr << "compiling answer..." << endl;
		{
			stringstream ss;
			ss << "cd " << testcase_dir << " && g++ answer.cpp -std=c++11 -O2 -o answer.exe";
			string cmd = ss.str();
			if(system(cmd.c_str()) != 0) {
				cerr << "an error occurred" << endl;
				abort();
			}
		}
	}

	cerr << "deleting files..." << endl;
	{
		string cmd = string("del \"") + testcase_dir + "\\*" + testcase_postfix + "\"";
		if(system(cmd.c_str()) != 0) {
			cerr << "an error occurred" << endl;
			abort();
		}
	}
	int numcases = (int)testcases.size();
	cerr << "number of test cases: " << numcases << endl;
	set<string> testnames;
	for(int ii = 0; ii < numcases; ii ++) {
		TestCase &t = *testcases[ii];

		const string testname = t.getName();
		cerr << "generating " << testname << "..." << endl;

		if(testnames.count(testname)) {
			cerr << "testname duplication: " << testname << endl;
			abort();
		}
		testnames.count(testname);

		//テストケース生成で使う乱数の初期化
		RNG rng;
		rep(i, problem_name.size())
			rng.addSeed(problem_name[i]);
		rep(i, testname.size())
			rng.addSeed(testname[i]);
		rep(i, 100) rng();

		//テストケース生成
		t.generate(rng);

		//ファイルへ書き込む
		string filename = string(testcase_dir) + "\\" + testname + testcase_postfix;
		fstream f(filename.c_str(), ios_base::out);
		if(!f) {
			cerr << "can't open the file" << endl;
			abort();
		}

		t.write(f);

		f.close();

		if(1) {
			cerr << "run " << testname << "..." << endl;
			string run_cmd = string(testcase_dir) + "\\answer.exe < \"" + filename + "\"";
			TIMEIT("run total") {
				if(system(run_cmd.c_str()) != 0) {
					cerr << "an error occurred" << endl;
					abort();
				}
			}
		}

		cerr << endl;
	}


	return 0;
}

#endif