#include <bits/stdc++.h>
using namespace std;
//#define cerr if (false) cerr
#define db(x) cerr << #x << "=" << x << endl
#define db2(x, y) cerr << #x << "=" << x << "," << #y << "=" << y << endl
#define db3(x, y, z) cerr << #x << "=" << x << "," << #y << "=" << y << "," << #z << "=" << z << endl
#define dbv(v) cerr << #v << "="; for (auto _x : v) cerr << _x << ", "; cerr << endl
#define dba(a, n) cerr << #a << "="; for (int _i = 0; _i < (n); ++_i) cerr << a[_i] << ", "; cerr << endl
template <typename A, typename B>
ostream& operator<<(ostream& os, const pair<A, B>& x) {
	return os << "(" << x.first << "," << x.second << ")";
}
typedef long long ll;
typedef long double ld;
const int MOD = 1000000007;
struct Mint {
	int val;
	Mint() { val = 0; }
	Mint(const ll& x) {
		val = (-MOD <= x && x < MOD) ? x : x % MOD;
		if (val < 0) val += MOD;
	}
	template <typename U>
	explicit operator U() const { return (U)val; }
	friend bool operator==(const Mint& a, const Mint& b) { return a.val == b.val; }
	friend bool operator!=(const Mint& a, const Mint& b) { return !(a == b); }
	friend bool operator<(const Mint& a, const Mint& b) { return a.val < b.val; }
	Mint& operator+=(const Mint& m) { if ((val += m.val) >= MOD) val -= MOD; return *this; }
	Mint& operator-=(const Mint& m) { if ((val -= m.val) < 0) val += MOD; return *this; }
	Mint& operator*=(const Mint& m) { val = (ll)val * m.val % MOD; return *this; }
	friend Mint modex(Mint a, ll p) {
		assert(p >= 0);
		Mint ans = 1;
		for (; p; p >>= 1, a *= a) if (p & 1) ans *= a;
		return ans;
	}
	Mint& operator/=(const Mint& m) { return *this *= modex(m, MOD - 2); }
	Mint& operator++() { return *this += 1; }
	Mint& operator--() { return *this -= 1; }
	Mint operator++(int) { Mint result(*this); *this += 1; return result; }
	Mint operator--(int) { Mint result(*this); *this -= 1; return result; }
	Mint operator-() const { return Mint(-val); }
	friend Mint operator+(Mint a, const Mint& b) { return a += b; }
	friend Mint operator-(Mint a, const Mint& b) { return a -= b; }
	friend Mint operator*(Mint a, const Mint& b) { return a *= b; }
	friend Mint operator/(Mint a, const Mint& b) { return a /= b; }
	friend ostream& operator<<(ostream& os, const Mint& x) { return os << x.val; }
};
const int MAXN = 5005;
Mint fac[MAXN], invfac[MAXN];
Mint binom(ll a, ll b) {
	if (b < 0) return 0;
	if (b > a) return 0;
	return fac[a] * invfac[b] * invfac[a - b];
}
vector<int> E[MAXN];
bool used[MAXN];
int sz;
bool dfs(int x) {
	if (x == -1) return true;
	if (used[x]) return false;
	++sz;
	used[x] = true;
	bool has = false;
	for (int y : E[x])
		if (dfs(y)) has = true;
	return has;
}
namespace fft {
struct cpx {
	double r, i;
	cpx(double r = 0, double i = 0) : r(r), i(i) {}
	cpx operator+(cpx b) { return {r + b.r, i + b.i}; }
	cpx operator-(cpx b) { return {r - b.r, i - b.i}; }
	cpx operator*(cpx b) { return {r * b.r - i * b.i, r * b.i + i * b.r}; }
	cpx operator/(double b) { return {r / b, i / b}; }
	inline cpx conj() { return {r, -i}; }
};
const double PI = acos(-1);
vector<int> rev(1, 0); // rev has [base] bits
vector<cpx> roots;
int base = 0;
void precompute(int nbase) {
	if (nbase <= base) return;
	rev.resize(1 << nbase);
	for (int i = 0; i < (1 << nbase); ++i) {
		rev[i] = rev[i >> 1] >> 1 | ((i & 1) << nbase - 1);
	}
	roots.resize(1 << nbase);
	for (; base < nbase; ++base) {
		int len = 1 << base;
		double angle = 2 * PI / (1 << (base + 1));
		for (int i = 0; i < len; ++i) {
			double cur = angle * i;
			roots[len + i] = cpx(cos(cur), sin(cur));
		}
	}
}
void fft(vector<cpx>& a, int n = -1) {
	if (n == -1) n = a.size();
	assert((n & (n - 1)) == 0);
	int zeros = __builtin_ctz(n);
	precompute(zeros);
	int shift = base - zeros;
	for (int i = 0; i < n; ++i) {
		if (i < (rev[i] >> shift)) swap(a[i], a[rev[i] >> shift]);
	}
	for (int len = 1; len < n; len <<= 1)
		for (int i = 0; i < n; i += 2 * len)
			for (int j = 0; j < len; ++j) {
				cpx u = a[i + j], v = a[i + j + len] * roots[len + j];
				a[i + j] = u + v;
				a[i + j + len] = u - v;
			}
}
vector<cpx> fa, fb;
vector<ll> multiply(const vector<int>& a, const vector<int>& b) {
	if (a.empty() || b.empty()) return {};
	int sz = a.size() + b.size() - 1;
	int n = sz == 1 ? 1 : 1 << (32 - __builtin_clz(sz - 1));
	if (n > fa.size()) fa.resize(n);
	for (int i = 0; i < n; ++i) {
		int x = i < a.size() ? a[i] : 0;
		int y = i < b.size() ? b[i] : 0;
		fa[i] = cpx(x, y);
	}
	fft(fa, n);
	cpx r(0, -0.5 / n);
	for (int i = 0; i <= (n >> 1); ++i) {
		int j = (n - i) & (n - 1);
		cpx x = (fa[i] + fa[j].conj()) / 2;
		cpx y = (fa[i] - fa[j].conj()) * r;
		fa[i] = x * y;
		if (i != j) fa[j] = fa[i].conj();
	}
	fft(fa, n);
	reverse(fa.begin() + 1, fa.begin() + n);
	vector<ll> res(sz);
	for (int i = 0; i < sz; ++i) {
		res[i] = llround(fa[i].r);
	}
	return res;
}
// Accepts Mint or int or ll in range [0..m-1].
template <typename T>
vector<T> multiply(const vector<T>& a, const vector<T>& b, int m, bool forceFast = false) {
	if (forceFast || m < (1 << 15)) {
		vector<int> aa(a.begin(), a.end()), bb(b.begin(), b.end());
		auto res = multiply(aa, bb);
		vector<T> ret(res.size());
		for (int i = 0; i < res.size(); ++i) ret[i] = res[i] % m;
		return ret;
	}
	if (a.empty() || b.empty()) return {};
	int eq = (a == b);
	int sz = a.size() + b.size() - 1;
	int n = sz == 1 ? 1 : 1 << (32 - __builtin_clz(sz - 1));
	if (n > fa.size()) fa.resize(n);
	for (int i = 0; i < a.size(); ++i) {
		fa[i] = cpx(int(a[i]) & ((1 << 15) - 1), int(a[i]) >> 15);
	}
	fill(fa.begin() + a.size(), fa.begin() + n, cpx(0, 0));
	fft(fa, n);
	if (n > fb.size()) fb.resize(n);
	if (eq) {
		copy(fa.begin(), fa.begin() + n, fb.begin());
	} else {
		for (int i = 0; i < b.size(); ++i) {
			fb[i] = cpx(int(b[i]) & ((1 << 15) - 1), int(b[i]) >> 15);
		}
		fill(fb.begin() + b.size(), fb.begin() + n, cpx(0, 0));
		fft(fb, n);
	}
	double ratio = 0.5 / n;
	cpx r2(0, -0.5);
	cpx r3(ratio, 0);
	cpx r4(0, -ratio);
	cpx r5(0, 1);
	for (int i = 0; i <= (n >> 1); ++i) {
		int j = (n - i) & (n - 1);
		cpx a1 = (fa[i] + fa[j].conj()) / 2;
		cpx a2 = (fa[i] - fa[j].conj()) * r2;
		cpx b1 = (fb[i] + fb[j].conj()) * r3;
		cpx b2 = (fb[i] - fb[j].conj()) * r4;
		cpx x = a1 * b1;
		cpx y = a2 * b2;
		fa[i] = x + y * r5;
		fb[i] = a1 * b2 + a2 * b1;
		if (i != j) {
			fa[j] = x.conj() + y.conj() * r5;
			fb[j] = fb[i].conj();
		}
	}
	fft(fa, n);
	fft(fb, n);
	reverse(fa.begin() + 1, fa.begin() + n);
	reverse(fb.begin() + 1, fb.begin() + n);
	vector<T> res(sz);
	for (int i = 0; i < sz; ++i) {
		ll aa = llround(fa[i].r);
		ll bb = llround(fb[i].r);
		ll cc = llround(fa[i].i);
		res[i] = (aa + ((bb % m) << 15) + ((cc % m) << 30)) % m;
	}
	return res;
}
/*
vector<Mint> multiply(const vector<Mint>& a, const vector<Mint>& b, bool forceFast = false) {
	return multiply(a, b, MOD, forceFast);
}
*/
}; // namespace fft
namespace polynomial {
// Leading zeros are not possible.
vector<Mint> subtract(vector<Mint> a, const vector<Mint>& b) {
	if (a.size() < b.size()) {
		a.resize(b.size());
	}
	for (int i = 0; i < b.size(); ++i) {
		a[i] -= b[i];
	}
	while (!a.empty() && a.back() == 0) a.pop_back();
	return a;
}
// Leading zeros are not possible.
vector<Mint> add(vector<Mint> a, const vector<Mint>& b) {
	if (a.size() < b.size()) {
		a.resize(b.size());
	}
	for (int i = 0; i < b.size(); ++i) {
		a[i] += b[i];
	}
	while (!a.empty() && a.back() == 0) a.pop_back();
	return a;
}
// Invert a power series mod x^n, where f[0] != 0
// P[2n] = P[n](2 - f * P[n]) mod x^(2n)
// O(n log n)
vector<Mint> inverse(const vector<Mint>& a) {
#ifdef NTT
	assert(!a.empty() && a[0] != 0);
	int n = a.size();
	vector<Mint> b = {1 / a[0]};
	for (int m = 2; b.size() < n; m <<= 1) {
		vector<Mint> aCut(a.begin(), a.begin() + min((int)a.size(), m));
		aCut.resize(2 * m);
		b.resize(2 * m);
		fft::fft(aCut);
		fft::fft(b);
		for (int i = 0; i < 2 * m; ++i) {
			b[i] *= 2 - aCut[i] * b[i];
		}
		fft::fft(b);
		reverse(b.begin() + 1, b.end());
		b.resize(m);
		Mint inv = Mint(1) / (2 * m);
		for (int i = 0; i < m; ++i) b[i] *= inv;
	}
	b.resize(n);
	return b;
#else
	assert(!a.empty() && a[0] != 0);
	int n = a.size();
	vector<Mint> b = {1 / a[0]};
	for (int m = 2; b.size() < n; m <<= 1) {
		vector<Mint> aCut(a.begin(), a.begin() + min((int)a.size(), m));
		auto x = fft::multiply(aCut, fft::multiply(b, b, MOD), MOD);
		b.resize(m);
		for (int i = 0; i < m; ++i) {
			b[i] += b[i] - x[i];
		}
	}
	b.resize(n);
	return b;
#endif
}
// Polynomial division: divide A (<= deg n) by B (deg m)
// A(x) = Q(x)B(x) + R(x), deg Q <= n-m. Substitute x = 1/x.
// x^nA(1/x) = (x^(n-m)Q(1/x))(x^mB(1/x)) + x^(n-m+1)(x^(m-1)R(1/x))
// Note x^deg * poly(1/x) = poly_reverse(x).
// So A_rev(x) = Q_rev(x)B_rev(x) mod x^(n-m+1).
// So we can recover Q, since deg Q_rev <= n-m.
// Leading zeros are not possible if A has no leading zeros.
// O(n log n)
vector<Mint> divide(vector<Mint> a, vector<Mint> b) {
	assert(!b.empty() && b.back() != 0);
	int n = a.size();
	int m = b.size();
	if (n < m) return {};
	reverse(a.begin(), a.end());
	reverse(b.begin(), b.end());
	b.resize(n - m + 1);
	a = fft::multiply(a, inverse(b), MOD);
	a.erase(a.begin() + n - m + 1, a.end());
	reverse(a.begin(), a.end());
	return a;
}
// O(n log n) where n = deg a. Leading zeros are not possible.
vector<Mint> remainder(const vector<Mint>& a, const vector<Mint>& b) {
	return subtract(a, fft::multiply(b, divide(a, b), MOD));
}
// Calculate integral with c=0. Leading zeros are possible!
vector<Mint> integral(const vector<Mint>& a) {
	vector<Mint> ret(a.size() + 1);
	for (int i = 1; i < ret.size(); ++i) {
		ret[i] = a[i - 1] / i;
	}
	return ret;
}
// Leading zeros are possible!
vector<Mint> derivative(const vector<Mint>& a) {
	if (a.empty()) return a;
	vector<Mint> ret(a.size() - 1);
	for (int i = 0; i + 1 < a.size(); ++i)
		ret[i] = a[i + 1] * (i + 1);
	return ret;
}
// Calculate log(f) mod x^n, where f[0] == 1. ret[0] == 0.
// O(n log n)
vector<Mint> logarithm(const vector<Mint>& a) {
	assert(!a.empty() && a[0] == 1);
	auto ret = integral(fft::multiply(derivative(a), inverse(a), MOD));
	ret.resize(a.size());
	return ret;
}
// Calculate exp(f) mod x^n, where f[0] == 0. ret[0] == 1.
// P[2n] = P[n](1 - log(P[n]) + f) mod x^(2n)
// O(n log n)
vector<Mint> exponent(const vector<Mint>& a) {
	assert(!a.empty() && a[0] == 0);
	int n = a.size();
	vector<Mint> b = {1};
	for (int m = 2; b.size() < n; m <<= 1) {
		vector<Mint> x(a.begin(), a.begin() + min((int)a.size(), m));
		++x[0];
		auto oldB = b;
		b.resize(m);
		b = fft::multiply(oldB, subtract(x, logarithm(b)), MOD);
		b.resize(m);
	}
	b.resize(n);
	return b;
}
// Multipoint evaluation: Evaluate A (deg d) at n points xs modulo MOD.
// Key idea: A(x_i) = A(x) mod (x - x_i).
// Calculate A mod (x-x_l)...(x-x_r) for range [l,r].
// O(n log^2 n + d log d)
vector<Mint> multipoint(const vector<Mint>& a, vector<Mint> xs) {
	int n = xs.size();
	vector<vector<Mint>> prods(2 * n - 1);
	function<void(int, int, int)> getProds = [&](int sn, int s, int e) {
		if (s == e) {
			prods[sn] = {-xs[s], 1};
			return;
		}
		int m = (s + e) >> 1, lsn = sn + 1, rsn = sn + ((m - s + 1) << 1);
		getProds(lsn, s, m);
		getProds(rsn, m + 1, e);
		prods[sn] = fft::multiply(prods[lsn], prods[rsn], MOD);
	};
	getProds(0, 0, n - 1);
	function<void(int, int, int, vector<Mint>)> calc = [&](int sn, int s, int e, vector<Mint> f) {
		if (s == e) {
			xs[s] = f.empty() ? 0 : f[0];
			return;
		}
		int m = (s + e) >> 1, lsn = sn + 1, rsn = sn + ((m - s + 1) << 1);
		calc(lsn, s, m, remainder(f, prods[lsn]));
		calc(rsn, m + 1, e, remainder(f, prods[rsn]));
	};
	calc(0, 0, n - 1, a);
	return xs;
}
// Given a polynomial A and generator g, evaluate Y[k] = A(g^k mod MOD) for all 0 <= k < n.
// A(g^(2k)) = sum_i a_i * g^(2ki)
// A(g^(2k+1)) = sum_i a_i * g^i * g^(2ki)
// Chirp z-transform: g^(2ki) = g^(kk) g^(ii) g^(-(k-i)^2); so we can convolve.
// O(MOD + (n + deg A) log (n + deg A))
vector<Mint> chirpz(vector<Mint> a, Mint g, int n) {
	int d = a.size();
	int phi = MOD - 1;
	vector<Mint> pw(phi + 1);
	pw[0] = 1;
	for (int i = 1; i <= phi; ++i) pw[i] = pw[i - 1] * g;
	auto chirpzEven = [&](const vector<Mint>& b, int m) {
		vector<Mint> u(d), v(m + d - 1);
		for (int i = 0; i < d; ++i) {
			u[i] = b[i] * pw[(ll)i * i % phi];
		}
		for (int i = -(d - 1); i < m; ++i) {
			v[d - 1 + i] = pw[phi - (ll)i * i % phi];
		}
		auto z = fft::multiply(u, v, MOD);
		vector<Mint> ret(m);
		for (int i = 0; i < m; ++i) {
			ret[i] = pw[(ll)i * i % phi] * z[d - 1 + i];
		}
		return ret;
	};
	auto even = chirpzEven(a, (n + 1) >> 1);
	for (int i = 0; i < d; ++i) a[i] *= pw[i % phi];
	auto odd = chirpzEven(a, n >> 1);
	vector<Mint> ret(n);
	for (int i = 0; i < n; ++i) {
		ret[i] = (i & 1) ? odd[i >> 1] : even[i >> 1];
	}
	return ret;
}
// Compute A^p mod B in O(b log b log p).
vector<Mint> powerMod(const vector<Mint>& a, ll p, const vector<Mint>& b) {
	if (!p) return {1};
	vector<Mint> c;
	if (p & 1) {
		c = fft::multiply(powerMod(a, p - 1, b), a, MOD);
	} else {
		c = powerMod(a, p >> 1, b);
		c = fft::multiply(c, c, MOD);
	}
	return remainder(c, b);
}
// Compute A^p mod x^n in O(n log n) using exp/log.
vector<Mint> power(vector<Mint> a, ll p, int n) {
	if (!p) return {1};
	assert(a.size() <= n);
	int m = 0;
	while (m < a.size() && a[m] == 0) ++m;
	if (m * p >= n) return {};
	a.erase(a.begin(), a.begin() + m);
	Mint inv = 1 / a[0];
	Mint c = modex(a[0], p);
	for (Mint& x : a) x *= inv;
	a = logarithm(a);
	Mint q = p;
	for (Mint& x : a) x *= q;
	a = exponent(a);
	a.insert(a.begin(), min(n - (ll)a.size(), m * p), 0);
	for (Mint& x : a) x *= c;
	return a;
}
vector<Mint> multiplyAll(int s, int e, const vector<vector<Mint>>& polys) {
	if (s == e) return polys[s];
	int m = (s + e) >> 1;
	auto a = multiplyAll(s, m, polys);
	auto b = multiplyAll(m + 1, e, polys);
	return fft::multiply(a, b, MOD);
}
}; // namespace polynomial
using namespace polynomial;
int main() {
	// Mint version.
	fac[0] = 1;
	for (int i = 1; i < MAXN; ++i)
		fac[i] = fac[i - 1] * i;
	invfac[MAXN - 1] = 1 / fac[MAXN - 1];
	for (int i = MAXN - 2; i >= 0; --i)
		invfac[i] = invfac[i + 1] * (i + 1);
	int n;
	scanf("%d", &n);
	int unknown = 0;
	for (int i = 1; i <= n; ++i) {
		int p;
		scanf("%d", &p);
		if (p == -1) ++unknown;
		E[i].push_back(p);
		if (p != -1) E[p].push_back(i);
	}
	vector<vector<Mint>> polys;
	Mint ans = 0;
	for (int i = 1; i <= n; ++i) {
		if (used[i]) continue;
		sz = 0;
		if (dfs(i)) {
			polys.push_back({1, sz});
			ans += (sz - 1) * modex(Mint(n - 1), unknown - 1);
		} else {
			ans += modex(Mint(n - 1), unknown);
		}
	}
	auto res = polys.empty() ? vector<Mint>{1} : multiplyAll(0, polys.size() - 1, polys);
	for (int k = 2; k <= min((int)res.size() - 1, unknown); ++k) {
		ans += res[k] * fac[k - 1] * modex(Mint(n - 1), unknown - k);
	}
	ans = n * modex(Mint(n - 1), unknown) - ans;
	printf("%d\n", ans.val);
}


Measured memory_bytes: 4988