#include <cstdio>
#include <cassert>
#include <ctime>

#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <utility>

#define rep(i, n) for (int i = 0; i < int(n); ++i)
#define rep3(i, m, n) for (int i = int(m); i < int(n); ++i)

using namespace std;
using uint64 = unsigned long long;
using uint32 = unsigned;

using Pair = pair<int, int>;

const int N_MAX = 4000;
const int MAT_SIZE = 4321;

const uint32 MOD = 1000000007;

Pair edges[N_MAX];

int vin[N_MAX];
int vout[N_MAX];
int conv[N_MAX];
bool deleted[N_MAX];

uint32 mat[MAT_SIZE][MAT_SIZE];
uint32 mat2[MAT_SIZE][MAT_SIZE];

uint32 facts[N_MAX + 1];

uint32 pow_mod(uint32 b, uint32 e) {
  uint32 ret = 1;
  while (e) {
    if (e & 1) ret = uint64(ret) * b % MOD;
    e >>= 1;
    b = uint64(b) * b % MOD;
  }
  return ret;
}

uint32 mod_inv(uint32 a) {
  return pow_mod(a, MOD - 2);
}

void init() {
  facts[0] = facts[1] = 1;
  for (uint32 j = 2; j <= N_MAX; ++j) {
    facts[j] = uint64(facts[j - 1]) * j % MOD;
  }
}

uint32 determinant(int n) {
  uint32 ret = 1;
  rep(i, n) {
    int j = i;
    if (mat[j][i] == 0) {
      for (++j ; j < n; ++j) if (mat[j][i] != 0) {
        rep3(k, i, n) swap(mat[i][k], mat[j][k]);
        break;
      }
      if (j == n) return 0;
      ret = MOD - ret;
    }
    ret = uint64(ret) * mat[i][i] % MOD;
    uint32 inv = mod_inv(MOD - mat[i][i]);
    rep3(j, i + 1, n) {
      uint64 c = mat[j][i];
      if (c == 0) continue;
      c = c * inv % MOD;
      rep3(k, i, n) mat[j][k] = (mat[j][k] + c * mat[i][k]) % MOD;
    }
  }
  return ret;
}

int calc(int N, int M, int vs, int ve) {
  rep(i, N) rep(j, N) mat[i][j] = (i == j ? N - 1 : MOD - 1);
  rep(i, M) {
    int vf = edges[i].first;
    int vt = edges[i].second;
    if (vf == vt) continue;
    mat[vf][vt] = 0;
    mat[vt][vt] -= 1;
  }
  if (vs >= 0) {
    mat[ve][vs] = (mat[ve][vs] + MOD - 1) % MOD;
    mat[vs][vs] += 1;
    vout[ve] += 1;
    vin[vs] += 1;
  }

  // semi-naive
  int v_removed = N - 1;
  int m_count = 0;
  rep3(i, 1, N) {
    if (mat[i][0] == 0) continue;
    int s = 0;
    rep(j, N) {
      s += mat[i][j] == MOD - 1;
    }
    if (s > m_count) {
      m_count = s;
      v_removed = i;
    }
  }
  rep(i, N) mat[0][i] = (MOD - mat[v_removed][i]) % MOD;
  if (v_removed < N - 1) {
    rep(i, N) rep3(j, v_removed, N - 1) mat[i][j] = mat[i][j + 1];
    rep3(i, v_removed, N - 1) rep(j, N) mat[i][j] = mat[i + 1][j];
  }
  rep3(i, 1, N - 1) rep(j, N - 1) {
    mat[i][j] = (mat[i][j] + mat[0][j]) % MOD;
  }
  // uint32 non_zero_count = 0;
  // rep(i, N - 1) rep(j, N - 1) non_zero_count += mat[i][j] != 0;
  // printf("[%u %u]\n", non_zero_count, v_removed);

  uint32 ret = determinant(N - 1);
  rep(i, N) {
    assert(vout[i] >= 1);
    ret = uint64(ret) * facts[vout[i] - 1] % MOD;
  }
  if (vs == -1) {
    ret = uint64(ret) * (N * N - M) % MOD;
  }
  return ret;
}

int solve(int N, int M) {
  rep(i, N) vin[i] = vout[i] = N, deleted[i] = false;

  rep(i, M) {
    int f = edges[i].first;
    int t = edges[i].second;
    vout[f] -= 1;
    vin[t] -= 1;
  }
  
  int v_total = 0;
  int vs = -1, ve = -1;
  rep(v, N) {
    if (vin[v] == vout[v]) {
      if (vin[v] == 0) {
        deleted[v] = true;
        continue;
      }
    } else if (vin[v] == vout[v] + 1) {
      if (ve != -1) return 0;
      ve = v_total;
    } else if (vin[v] + 1 == vout[v]) {
      if (vs != -1) return 0;
      vs = v_total;
    } else {
      return 0;
    }
    vin[v_total] = vin[v];
    vout[v_total] = vout[v];
    conv[v] = v_total++;
  }

  if (v_total == 0) return 1;

  int e_total = 0;
  rep(e, M) {
    int f = edges[e].first;;
    int t = edges[e].second;
    if (deleted[f] || deleted[t]) continue;
    edges[e_total++] = Pair(conv[f], conv[t]);
  }

  /* shuffle test */
  rep(i, v_total) conv[i] = i;
  srand(time(nullptr));
  for (int i = v_total - 1; i > 0; --i) {
    int j = rand() % (i + 1);
    swap(conv[i], conv[j]);
    swap(vout[conv[i]], vout[conv[j]]);
    swap(vin[conv[i]], vin[conv[j]]);
  }
  if (vs != -1) {
    vs = conv[vs];
    ve = conv[ve];
  }
  rep(e, e_total) edges[e] = Pair(conv[edges[e].first], conv[edges[e].second]);

  return calc(v_total, e_total, vs, ve);
}

int main() {
  init();

  int N, M; 
  while (~scanf("%d %d", &N, &M)) {
    assert(N <= N_MAX && M <= N_MAX && M <= N * N);

    int a, b;
    rep(i, M) scanf("%d %d", &a, &b), edges[i] = Pair(a - 1, b - 1);

    int ans = solve(N, M);
    printf("%d\n", ans);
  }
  return 0;
}