#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef complex<double> P;
typedef pair<int,int> pii;
#define REP(i,n) for(ll i=0;i<n;++i)
#define REPR(i,n) for(ll i=1;i<n;++i)
#define FOR(i,a,b) for(ll i=a;i<b;++i)

#define DEBUG(x) cout<<#x<<": "<<x<<endl
#define DEBUG_VEC(v) cout<<#v<<":";REP(i,v.size())cout<<" "<<v[i];cout<<endl
#define ALL(a) (a).begin(),(a).end()

#define MOD (ll)(1e9+7)
#define ADD(a,b) a=((a)+(b))%MOD
#define FIX(a) ((a)%MOD+MOD)%MOD

#define V_MAX 4010

struct UF{
  vi data;
  UF(int size):data(size,-1){}
  int root(int a){
    return data[a]<0 ? a : data[a]=root(data[a]);
  }
  void unite(int a,int b){
    a=root(a);
    b=root(b);
    if(a!=b){
      if(data[b]<data[a])swap(a,b);
      data[a] += data[b];
      data[b] = a;
    }
  }
  bool same(int a,int b){
    return root(a) == root(b);
  }
  int size(int a){
    return -data[root(a)];
  }
};

ll fact[V_MAX];
void fact_init(){
  fact[0] = 1;
  REPR(i,V_MAX){
    fact[i] = i*fact[i-1];
    fact[i] %= MOD;
  }
}

int inv(ll a){
  ll t = MOD-2;
  ll res = 1;
  while(t){
    if(t&1==1){
      res = res*a%MOD;
    }
    t >>= 1;
    a = a*a%MOD;
  }
  return res;
}

ll determinant(int dim,vector<vl> mat){
  if(dim==0)return 1;
  REPR(i,dim)REP(j,dim){
    mat[i][j] -= mat[0][j];
  }
  REP(i,dim)REP(j,dim)mat[i][j] = FIX(mat[i][j]);
  ll result = 1;
  REP(i,dim-1){
    bool flag = false;
    FOR(j,i,dim){
      if(mat[j][i]!=0){
        if(i!=j){
          swap(mat[j],mat[i]);
          result = (-result + MOD)%MOD;
        }
        flag = true;
        break;
      }
    }
    if(!flag)return 0;
    ll iv = inv(mat[i][i]);
    FOR(j,i+1,dim){
      if(mat[j][i]==0)continue;
      ll rate = mat[j][i] * iv;
      rate %= MOD;
      FOR(k,i,dim){
        mat[j][k] -= (mat[i][k]*rate)%MOD;
        mat[j][k] = FIX(mat[j][k]);
      }
    }
    result *= mat[i][i];
    result %= MOD;
  }
  result *= mat[dim-1][dim-1];
  result %= MOD;
  return result;
}

ll solve(){
  int n,m;
  cin >> n >> m;
  vi indeg(n,n),outdeg(n,n);
  pii edges[m];
  REP(i,m){
    int a,b;
    cin >> a >> b;
    --a;--b;
    edges[i] = make_pair(a,b);
    outdeg[a]--;
    indeg[b]--;
  }
  // no edge
  if(n*n==m)return 1;

  // check eulerian
  int beg=-1,end=-1;
  REP(i,n){
    if(indeg[i]==outdeg[i])continue;
    if(beg==-1 && indeg[i]-1==outdeg[i]){
      beg=i;
      continue;
    }
    if(end==-1 && indeg[i]==outdeg[i]-1){
      end=i;
      continue;
    }
    return 0;
  }
  if(beg!=-1 && end==-1)return 0;
  if(beg==-1 && end!=-1)return 0;
  if(beg!=-1 && end!=-1){
    outdeg[beg]++;
    indeg[end]++;
  }

  // check & delete isolation
  sort(edges,edges+m);

  UF uf = UF(n);
  int root = -1;
  int iter = 0;
  REP(i,n){
    if(root==-1 && indeg[i]>0) root = i;
    REP(j,n){
      if(iter<m && edges[iter].first == i && edges[iter].second == j){
        ++iter;
      }else{
        uf.unite(i,j);
      }
    }
  }
  if(root==-1)return 1;
  int mp[n];
  int t = 0;
  vi deg(n);
  REP(i,n){
    if(!uf.same(root,i)){
      if(indeg[i]>0){
        return 0;
      }
      mp[i] = -1;
    }else{
      mp[i] = t;
      deg[t] = outdeg[i];
      ++t;
    }
  }
  REP(i,m){
    edges[i].first = mp[edges[i].first];
    edges[i].second = mp[edges[i].second];
  }
  vector<vl> Q(t,vl(t,-1));
  REP(i,t){
    Q[i][i] += deg[i];
  }
  REP(i,m){
    if(edges[i].first != -1 && edges[i].second != -1)
      Q[edges[i].first][edges[i].second] += 1;
  }
  if(beg!=-1 && end!=-1){
    beg = mp[beg];
    end = mp[end];
    if(beg==-1 || end==-1)return 0;
    Q[end][beg] -= 1;
  }

  // BEST theorem
  // ec(G) = t_w(G) PI_{v in V}(deg(v)-1)!
  // matrix tree theorem
  // t_w(G) = (determinant of Laplacian matrix)
  fact_init();
  ll result = 1;
  REP(i,t){
    result *= fact[deg[i]-1];
    result %= MOD;
  }
  result *= determinant(t-1,Q);
  result %= MOD;
  if(beg==-1 && end==-1){
    result *= n*n-m;
    result %= MOD;
  }
  return result;
}

int main(){
  // counting eulerian circuit
  cout << solve() << endl;
  return 0;
}