#include <bits/stdc++.h>
#define GET_MACRO(a, b, c, NAME, ...) NAME
#define rep(...) GET_MACRO(__VA_ARGS__, rep3, rep2)(__VA_ARGS__)
#define rep2(i, a) rep3 (i, 0, a)
#define rep3(i, a, b) for (int i = (a); i < (b); i++)
#define repr(...) GET_MACRO(__VA_ARGS__, repr3, repr2)(__VA_ARGS__)
#define repr2(i, a) repr3 (i, 0, a)
#define repr3(i, a, b) for (int i = (b) - 1; i >= (a); i--)
template<class T1, class T2> inline bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> inline bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
using namespace std;
typedef long long ll;

const ll mod = 1e9 + 7;

struct SparseDeterminant {
	unsigned long long xor64() {
		static unsigned long long x = time(NULL);
		x ^= x << 13; x ^= x >> 7; x ^= x << 17;
		return x;
	}

	vector<ll> random_vector(int n) {
		vector<ll> res(n);
		rep (i, n) res[i] = xor64() % mod;
		return res;
	}

	// Berlekamp-Massey algorithm
	vector<ll> minpoly(vector<ll> a) {
		const int N = a.size();
		vector<ll> b(N), c(N), t(N);
		b[0] = 1;
		c[0] = 1;
		int l = 0;
		int m = -1;
		for (int n = 0; n < N; n++) {
			int d = 0;
			for (int i = 0; i <= l; i++) {
				(d += c[i] * a[n - i]) %= mod;
			}
			if (d == 1) {
				t = c;
				int N_M = n - m;
				for (int j = 0; j < N - N_M; j++) {
					(c[N_M + j] += b[j]) %= mod;
				}
				if (l <= n / 2) {
					l = n + 1 - l;
					m = n;
					b = t;
				}
			}
		}
		return c;
	}

	vector<ll> minpoly_vector(vector<vector<ll>> b) {
		// TODO: implement
	}

	vector<ll> minpoly_matrix(vector<vector<ll>> A) {
		// TODO: implement
	}
};


ll F[5050];

ll modpow(ll a, ll b, ll mod) {
	ll res = 1;
	while (b) {
		if (b & 1) (res *= a) %= mod;
		(a *= a) %= mod;
		b /= 2;
	}
	return res;
}

ll modinv(ll a, ll mod) {
	return modpow(a, mod - 2, mod);
}

void show(vector<vector<ll>> A) {
	rep (i, A.size()) {
		rep (j, A[0].size()) {
			cout << A[i][j] << " ";
		}
		cout << endl;
	}
	cout << endl;
}

ll det(vector<vector<ll>> A) {
	int n = A.size();
	rep (j, n) {
		rep (i, j, n) {
			if (A[i][j] != 0) {
				swap(A[j], A[i]);
				break;
			}
		}
		ll inv = modinv(A[j][j], mod);
		rep (i, j + 1, n) {
			repr (k, j, n) {
				A[i][k] -= A[i][j] * A[j][k] % mod * inv % mod;
				if (A[i][k] < 0) A[i][k] += mod;
			}
		}
	}
	ll res = 1;
	rep (i, n) (res *= A[i][i]) %= mod;
	return res;
}

vector<vector<ll>> remove(vector<vector<ll>> A, int y, int x) {
	int N = A.size();
	rep (i, N) {
		rep (j, x, N - 1) A[i][j] = A[i][j + 1];
	}
	rep (j, N) {
		rep (i, y, N - 1) A[i][j] = A[i + 1][j];
	}
	rep (i, N) A[i].resize(N - 1);
	A.resize(N - 1);
	return A;
}


int main() {
	F[0] = 1;
	rep (i, 1, 5050) F[i] = i * F[i - 1] % mod;

	int N, M;
	cin >> N >> M;

	vector<vector<ll>> L(N, vector<ll>(N)); // laplacian matrix
	vector<vector<ll>> A(N, vector<ll>(N)); // adjacent matrix
	vector<vector<ll>> D(N, vector<ll>(N)); // degree matrix
	vector<int> outdeg(N, N), indeg(N, N);

	rep (i, N) rep (j, N) A[i][j] = 1;
	rep (i, N) D[i][i] = N;

	rep (i, M) {
		int u, v;
		cin >> u >> v;
		u--; v--;
		D[u][u]--;
		A[u][v] = 0;
		outdeg[u]--;
		indeg[v]--;
	}

	int u = -1, v = -1;
	rep (i, N) {
		if (indeg[i] == outdeg[i] - 1) {
			if (u == -1) {
				u = i;
			} else {
				cout << 0 << endl;
				return 0;
			}
		} else if (indeg[i] == outdeg[i] + 1) {
			if (v == -1) {
				v = i;
			} else {
				cout << 0 << endl;
				return 0;
			}
		} else if (indeg[i] != outdeg[i]) {
			cout << 0 << endl;
			return 0;
		}
	}

	bool cycle = false;
	if (u == -1 && v == -1) {
		u = 0, v = 0;
		cycle = true;
	} else if (u != -1 ^ v != -1) {
		cout << 0 << endl;
		return 0;
	}

	rep (i, N) rep (j, N) L[i][j] = D[i][j] - A[i][j];

	auto Lw = remove(L, u, u);
	ll d = det(Lw);

	ll ans = d;

	rep (i, N) (ans *= F[outdeg[i] - 1]) %= mod;
	if (cycle) (ans *= N * N - M) %= mod;

	cout << ans << endl;

	return 0;
}