#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
#define rep(i,s,e) for(i64 (i) = (s);(i) < (e);(i)++)
#define all(x) x.begin(),x.end()

template<class T>
static inline std::vector<T> ndvec(size_t&& n, T val) { return std::vector<T>(n, std::forward<T>(val)); }
template<class... Tail>
static inline auto ndvec(size_t&& n, Tail&&... tail) {
  return std::vector<decltype(ndvec(std::forward<Tail>(tail)...))>(n, ndvec(std::forward<Tail>(tail)...));
}
template<class Cond> struct chain {
  Cond cond; chain(Cond cond) : cond(cond) {}
  template<class T> bool operator()(T& a, const T& b) const { if(cond(a, b)) { a = b; return true; } return false; }
};
template<class Cond> chain<Cond> make_chain(Cond cond) { return chain<Cond>(cond); }

#include <iostream>
using i64 = long long;
template<i64 M> struct modint { i64 a;
  constexpr modint(const i64 x = 0): a((x%M+M)%M){}
  constexpr i64 value() const { return a; }
  constexpr modint inv() const { return this->pow(M-2); }
  constexpr modint pow(i64 r) const {
    modint ans(1); modint aa = *this;
    while(r) { if(r & 1) ans *= aa; aa *= aa; r >>= 1; }
    return ans;
  }
  constexpr modint& operator=(const i64 r) { a = (r % M + M) % M; return *this; }
  constexpr modint& operator+=(const modint r) { a += r.a; if(a >= M) a -= M; return *this; }
  constexpr modint& operator-=(const modint r) { a -= r.a; if(a < 0) a += M; return *this; }
  constexpr modint& operator*=(const modint r) { a = a * r.a % M; return *this; }
  constexpr modint& operator/=(const modint r) { (*this) *= r.inv(); return *this; }
  constexpr modint operator+(const modint r) const { return modint(*this) += r; }
  constexpr modint operator-(const modint r) const { return modint(*this) -= r; }
  constexpr modint operator*(const modint r) const { return modint(*this) *= r; }
  constexpr modint operator/(const modint r) const { return modint(*this) /= r; }
  constexpr bool operator!=(const modint r) const { return this->value() != r.value(); }
};

template<const i64 M> std::ostream& operator<<(std::ostream& os, const modint<M>& m) { os << m.value(); return os; }

#include <vector>
#include <tuple>

struct union_find {
  std::vector<int> par;
  union_find(int N): par(N, -1) {}
  int root(int x) {
    return par[x] < 0 ? x : par[x] = root(par[x]);
  }
  std::tuple<int, int> unite(int x, int y) {
    x = root(x);
    y = root(y);
    if(x == y) return { -1, -1 };
    if(par[x] > par[y]) std::swap(x, y);
    par[x] += par[y];
    par[y] = x;
    return { x, y };
  }
  int size(int x) {
    return -par[root(x)];
  }
};

using fp = modint<(i64)(1e9 + 7)>;

#include <vector>
#include <tuple>
using i64 = long long;

template<class T>
std::tuple<std::vector<T>, std::vector<T>, std::vector<T>> build_factorial(int N) {
  std::vector<T> fact(N);
  std::vector<T> finv(N);
  std::vector<T> inv(N);
  fact[0] = T(1);
  for(int i = 1;i < N;i++) {
    fact[i] = fact[i - 1] * T(i);
  }
  finv[N - 1] = T(1) / fact[N - 1];
  for(int i = N - 1; i --> 0;) {
    finv[i] = finv[i + 1] * T(i + 1);
  }
  for(int i = 0;i < N;i++) {
    inv[i] = fact[i - 1] * finv[i];
  }
  return std::make_tuple(std::move(fact), std::move(finv), std::move(inv));
}

int main() {
  i64 N;
  cin >> N;

  auto [fact, finv, inv] = build_factorial<fp>(N + 1);
  auto binom = [&](i64 n, i64 k) {
    if(0 <= k && k <= n) return fact[n] * finv[k] * finv[n - k];
    return fp(0);
  };

  vector<i64> P(N);
  rep(i,0,N) {
    cin >> P[i];
    P[i]--;
  }
  i64 K = 0;
  union_find uf(N);
  rep(i,0,N) {
    if(P[i] != -2) {
      uf.unite(i, P[i]);
    }
    else {
      K++;
    }
  }

  set<i64> st;
  rep(i,0,N) {
    i64 r = uf.root(i);
    if(st.count(r)) continue;
    st.insert(r);
  }
  vector<i64> sz;
  rep(i,0,N) {
    if(P[i] == -2) {
      sz.push_back(uf.size(i));
    }
  }

  auto fdp = ndvec(N + 1, N + 1, fp(0));
  fdp[0][0] = 1;

  fp res(0);
  res += fp(st.size() - sz.size()) * fp(N - 1).pow(sz.size());

  rep(i,0,sz.size()) {
    res += fp(sz[i] - 1) * fp(N - 1).pow(sz.size() - 1);
    rep(j,0,N + 1) {
      if(j > 0) {
        fdp[i + 1][j] += fdp[i][j - 1] * sz[i];
      }
      fdp[i + 1][j] += fdp[i][j];
      if(j > 0 && j + 1 <= sz.size()) {
        res += fdp[i][j] * sz[i] * fact[j] * fp(N - 1).pow(sz.size() - j - 1);
      }
    }
  }
  cout << fp(N) * fp(N - 1).pow(sz.size()) - res << endl;
}


Measured memory_bytes: 199312