#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 C[V_MAX+1],B[V_MAX],T[V_MAX];
ll _N;

ll berlekamp_massey(ll s[]){
  ll N = 2*_N;
  ll L=0,m=1,b=1,Bsz=1;
  C[0]=1;
  B[0]=1;
  REP(n,N){
    ll d = s[n];
    REPR(i,L+1) d += C[i]*s[n-i]%MOD;
    d %= MOD;
    if(d==0){
      ++m;
    }else if(2*L <= n){
      // copy C to T
      // calc C from C&B
      // copy T to B
      // C->T, C<-C&B, B<-T
      REP(i,L+1)T[i]=C[i];
      ll rate = (MOD-(d*inv(b))%MOD)%MOD;
      REP(i,Bsz) C[m+i] = (C[m+i]+rate*B[i])%MOD;
      Bsz = L+1;
      L = n+1-L;
      REP(i,Bsz)B[i]=T[i];
      b = d;
      m = 1;
    }else{
      ll rate = (MOD-(d*inv(b))%MOD)%MOD;
      REP(i,Bsz) C[m+i] = (C[m+i]+rate*B[i])%MOD;
      ++m;
    }
  }
  return L+1;
}

typedef unsigned long ul;
ul xorshift(){
  static ul x = 123456789,
            y = 362436069,
            z = 521288629,
            w = 88675123;
  ul t;
  t = x^(x<<11);
  x = y;
  y = z;
  z = w;
  return w = (w^(w>>19))^(t^(t>>8));
}

pii edges[V_MAX+1];

// 対角成分はなんか値が入ってる
// それ以外の「ほとんど」は-1
// -1じゃない(0)のはたかだかM個
ll det(ll n,ll dmat[],ll m,ll beg,ll end){
  if(n<=0)return 1;
  _N = n;
  vl b(n),u(n);
  vl D(n); // diagonal matrix
  REP(i,n)while(b[i]==0)b[i] = xorshift()%MOD;
  REP(i,n)while(u[i]==0)u[i] = xorshift()%MOD;
  REP(i,n)while(D[i]==0)D[i] = xorshift()%MOD;
  ll a[2*n];
  a[0]=0;
  REP(j,n)a[0]+=u[j]*b[j]%MOD;
  a[0]%=MOD;
  REPR(i,2*n){
    // a_i = u^T * ( mat*D * (mat*D)^(i-1) * b)
    // 右から
    vl mb(n);
    ll sum = 0;
    REP(j,n){
      b[j] = (b[j]*D[j])%MOD;
      sum = (sum+b[j])%MOD;
    }
    REP(j,n){
      mb[j] = FIX(b[j]*dmat[j]-sum);
    }
    if(beg!=-1 && end!=-1 && beg<n && end<n){
      mb[end] = FIX(mb[end]-b[beg]);
    }
    REP(j,m){
      pii e = edges[j];
      if(e.first!=-1 && e.second!=-1 && e.first<n && e.second<n){
        mb[e.first] = FIX(mb[e.first]+b[e.second]);
      }
    }
    b = mb;
    a[i] = 0;
    REP(j,n){
      a[i] = FIX(a[i]+u[j]*b[j]);
    }
  }

  ll minimal = berlekamp_massey(a);
  if(minimal<n+1)return -1;
  ll ret = 1;
  ll detd = 1;
  REP(i,n)detd = (detd*D[i])%MOD;
  ret = C[n] * inv(detd) % MOD;
  if(n%2==1)ret = MOD-ret;
  ret %= MOD;
  return ret;
}

ll solve(){
  xorshift();
  int n,m;
  cin >> n >> m;
  vi indeg(n,n),outdeg(n,n);
  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]==outdeg[i]-1){
      beg=i;
      continue;
    }
    if(end==-1 && indeg[i]-1==outdeg[i]){
      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[end]++;
    indeg[beg]++;
  }

  // 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;
  ll 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];
  }
  if(beg!=-1 && end!=-1){
    beg = mp[beg];
    end = mp[end];
    if(beg==-1 || end==-1)return 0;
  }

  // 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;
  }
  ll dt = det(t-1,deg,m,beg,end);
  result *= dt;
  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;
}