#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;
}

vl berlekamp_massey(vl s){
  ll N = s.size();
  ll L=0,m=1,b=1,Bsz=1;
  vl C(N+1),B(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){
      vl T;REP(i,L+1)T.push_back(C[i]);
      ll rate = (MOD-(d*inv(b))%MOD)%MOD;
      REP(i,Bsz) C[m+i] = (C[m+i]+rate*B[i])%MOD;
      L = n+1-L;
      B = T;
      Bsz = B.size();
      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;
    }
  }
  C.resize(L+1);
  return C;
}

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

// 対角成分はなんか値が入ってる
// それ以外の「ほとんど」は-1
// -1じゃない(0)のはたかだかM個
ll det(ll n,vl dmat,vector<pii> edges,ll beg,ll end){
  if(n<=0)return 1;
  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;
  vl 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,edges.size()){
      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]);
    }
  }

  vl minimal = berlekamp_massey(a);
  if(a.size()<n+1)return det(n,dmat,edges,beg,end);
  ll ret = 1;
  ll detd = 1;
  REP(i,n)detd = (detd*D[i])%MOD;
  ret = minimal[n] * inv(detd) % MOD;
  if(n%2==1)ret = MOD-ret;
  ret %= MOD;
  return ret;
}

// 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(){
  xorshift();
  int n,m;
  cin >> n >> m;
  vi indeg(n,n),outdeg(n,n);
  vector<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]==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(ALL(edges));

  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;
  vl 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;
  }

  // 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){
  //   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;
  }
  ll dt = det(t-1,deg,edges,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;
}