//header
#ifdef LOCAL
    #include "cxx-prettyprint-master/prettyprint.hpp"
    #define debug(x) cout << x << endl
#else
    #define debug(...) 42
#endif
    #pragma GCC optimize("Ofast")
    #include <bits/stdc++.h>
    //types
    using namespace std;
    using ll = long long;
    using ul = unsigned long long;
    using ld = long double;
    typedef pair < ll , ll > Pl;        
    typedef pair < int, int > Pi;
    typedef vector<ll> vl;
    typedef vector<int> vi;
    template< typename T >
    using mat = vector< vector< T > >;
    template< int mod >
    struct modint {
        int x;

        modint() : x(0) {}

        modint(int64_t y) : x(y >= 0 ? y % mod : (mod - (-y) % mod) % mod) {}

        modint &operator+=(const modint &p) {
            if((x += p.x) >= mod) x -= mod;
            return *this;
        }

        modint &operator-=(const modint &p) {
            if((x += mod - p.x) >= mod) x -= mod;
            return *this;
        }

        modint &operator*=(const modint &p) {
            x = (int) (1LL * x * p.x % mod);
            return *this;
        }

        modint &operator/=(const modint &p) {
            *this *= p.inverse();
            return *this;
        }

        modint operator-() const { return modint(-x); }

        modint operator+(const modint &p) const { return modint(*this) += p; }

        modint operator-(const modint &p) const { return modint(*this) -= p; }

        modint operator*(const modint &p) const { return modint(*this) *= p; }

        modint operator/(const modint &p) const { return modint(*this) /= p; }

        bool operator==(const modint &p) const { return x == p.x; }

        bool operator!=(const modint &p) const { return x != p.x; }

        modint inverse() const {
            int a = x, b = mod, u = 1, v = 0, t;
            while(b > 0) {
            t = a / b;
            swap(a -= t * b, b);
            swap(u -= t * v, v);
            }
            return modint(u);
        }

        modint pow(int64_t n) const {
            modint ret(1), mul(x);
            while(n > 0) {
            if(n & 1) ret *= mul;
            mul *= mul;
            n >>= 1;
            }
            return ret;
        }

        friend ostream &operator<<(ostream &os, const modint &p) {
            return os << p.x;
        }

        friend istream &operator>>(istream &is, modint &a) {
            int64_t t;
            is >> t;
            a = modint< mod >(t);
            return (is);
        }

        static int get_mod() { return mod; }
    };
    //abreviations
    #define all(x) (x).begin(), (x).end()
    #define rall(x) (x).rbegin(), (x).rend()
    #define rep_(i, a_, b_, a, b, ...) for (int i = (a), max_i = (b); i < max_i; i++)
    #define rep(i, ...) rep_(i, __VA_ARGS__, __VA_ARGS__, 0, __VA_ARGS__)
    #define rrep_(i, a_, b_, a, b, ...) for (int i = (b-1), min_i = (a); i >= min_i; i--)
    #define rrep(i, ...) rrep_(i, __VA_ARGS__, __VA_ARGS__, 0, __VA_ARGS__)
    #define SZ(x) ((int)(x).size())
    #define pb(x) push_back(x)
    #define eb(x) emplace_back(x)
    #define mp make_pair
    #define print(x) cout << x << endl
    #define vsum(x) accumulate(x, 0LL)
    #define vmax(a) *max_element(all(a))
    #define vmin(a) *min_element(all(a))
    //functions
    ll gcd(ll a, ll b) { return b ? gcd(b, a%b) : a; }
    ll lcm(ll a, ll b) { return a/gcd(a, b)*b;}
    template<class T>bool chmax(T &a, const T &b) { if (a<b) { a=b; return 1; } return 0; }
    template<class T>bool chmin(T &a, const T &b) { if (b<a) { a=b; return 1; } return 0; }
    template< typename T >
    T mypow(T x, ll n) {
        T ret = 1;
        while(n > 0) {
            if(n & 1) (ret *= x);
            (x *= x);
            n >>= 1;
        }
        return ret;
    }
    ll modpow(ll x, ll n, const ll mod) {
        ll ret = 1;
        while(n > 0) {
            if(n & 1) (ret *= x);
            (x *= x);
            n >>= 1;
            x%=mod;
            ret%=mod;
        }
        return ret;
    }
    uint64_t my_rand(void) {
        static uint64_t x = 88172645463325252ULL;
        x = x ^ (x << 13); x = x ^ (x >> 7);
        return x = x ^ (x << 17);
    }
    //graph template
    template< typename T >
    struct edge {
        int src, to;
        T cost;

        edge(int to, T cost) : src(-1), to(to), cost(cost) {}

        edge(int src, int to, T cost) : src(src), to(to), cost(cost) {}

        edge &operator=(const int &x) {
            to = x;
            return *this;
        }
        operator int() const { return to; }
    };
    template< typename T >
    using Edges = vector< edge< T > >;
    template< typename T >
    using WeightedGraph = vector< Edges< T > >;
    using UnWeightedGraph = vector< vector< int > >;

//constant
#define INF 4001002003004005006LL
#define inf 1000000005
#define mod 1000000007LL
#define endl '\n'
typedef modint<mod> mint;
const long double eps = 0.001;
const long double PI  = 3.141592653589793;
//library
struct UnionFind {
        vector< int > data;

        UnionFind(int sz) {
                data.assign(sz, -1);
        }

        bool unite(int x, int y) {
                x = find(x), y = find(y);
                if(x == y) return (false);
                if(data[x] > data[y]) swap(x, y);
                data[x] += data[y];
                data[y] = x;
                return (true);
        }

        int find(int k) {
                if(data[k] < 0) return (k);
                return (data[k] = find(data[k]));
        }

        int size(int k) {
                return (-data[find(k)]);
        }
        bool same(int x, int y){
                return find(x)==find(y);
        }
};

template< typename T >
struct FormalPowerSeries : vector< T > {
    using vector< T >::vector;
    using P = FormalPowerSeries;

    using MULT = function< P(P, P) >;

    static MULT &get_mult() {
        static MULT mult = nullptr;
        return mult;
    }

    static void set_fft(MULT f) {
        get_mult() = f;
    }

    void shrink() {
        while(this->size() && this->back() == T(0)) this->pop_back();
    }

    P operator+(const P &r) const { return P(*this) += r; }

    P operator+(const T &v) const { return P(*this) += v; }

    P operator-(const P &r) const { return P(*this) -= r; }

    P operator-(const T &v) const { return P(*this) -= v; }

    P operator*(const P &r) const { return P(*this) *= r; }

    P operator*(const T &v) const { return P(*this) *= v; }

    P operator/(const P &r) const { return P(*this) /= r; }

    P operator%(const P &r) const { return P(*this) %= r; }

    P &operator+=(const P &r) {
        if(r.size() > this->size()) this->resize(r.size());
        for(int i = 0; i < r.size(); i++) (*this)[i] += r[i];
        return *this;
    }

    P &operator+=(const T &r) {
        if(this->empty()) this->resize(1);
        (*this)[0] += r;
        return *this;
    }

    P &operator-=(const P &r) {
        if(r.size() > this->size()) this->resize(r.size());
        for(int i = 0; i < r.size(); i++) (*this)[i] -= r[i];
        shrink();
        return *this;
    }

    P &operator-=(const T &r) {
        if(this->empty()) this->resize(1);
        (*this)[0] -= r;
        shrink();
        return *this;
    }

    P &operator*=(const T &v) {
        const int n = (int) this->size();
        for(int k = 0; k < n; k++) (*this)[k] *= v;
        return *this;
    }

    P &operator*=(const P &r) {
        if(this->empty() || r.empty()) {
            this->clear();
            return *this;
        }
        assert(get_mult() != nullptr);
        return *this = get_mult()(*this, r);
    }

    P &operator%=(const P &r) {
        return *this -= *this / r * r;
    }

    P operator-() const {
        P ret(this->size());
        for(int i = 0; i < this->size(); i++) ret[i] = -(*this)[i];
        return ret;
    }

    P &operator/=(const P &r) {
        if(this->size() < r.size()) {
            this->clear();
            return *this;
        }
        int n = this->size() - r.size() + 1;
        return *this = (rev().pre(n) * r.rev().inv(n)).pre(n).rev(n);
    }

    P pre(int sz) const {
        return P(begin(*this), begin(*this) + min((int) this->size(), sz));
    }

    P operator>>(int sz) const {
        if(this->size() <= sz) return {};
        P ret(*this);
        ret.erase(ret.begin(), ret.begin() + sz);
        return ret;
    }

    P operator<<(int sz) const {
        P ret(*this);
        ret.insert(ret.begin(), sz, T(0));
        return ret;
    }

    P rev(int deg = -1) const {
        P ret(*this);
        if(deg != -1) ret.resize(deg, T(0));
        reverse(begin(ret), end(ret));
        return ret;
    }

    P diff() const {
        const int n = (int) this->size();
        P ret(max(0, n - 1));
        for(int i = 1; i < n; i++) ret[i - 1] = (*this)[i] * T(i);
        return ret;
    }

    P integral() const {
        const int n = (int) this->size();
        P ret(n + 1);
        ret[0] = T(0);
        for(int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / T(i + 1);
        return ret;
    }

    // F(0) must not be 0
    P inv(int deg = -1) const {
        assert(((*this)[0]) != T(0));
        const int n = (int) this->size();
        if(deg == -1) deg = n;
        P ret({T(1) / (*this)[0]});
        for(int i = 1; i < deg; i <<= 1) {
            ret = (ret + ret - ret * ret * pre(i << 1)).pre(i << 1);
        }
        return ret.pre(deg);
    }

    // F(0) must be 1
    P log(int deg = -1) const {
        assert((*this)[0] == 1);
        const int n = (int) this->size();
        if(deg == -1) deg = n;
        return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
    }

    P sqrt(int deg = -1) const {
        const int n = (int) this->size();
        if(deg == -1) deg = n;

        if((*this)[0] == T(0)) {
            for(int i = 1; i < n; i++) {
                if((*this)[i] != T(0)) {
                    if(i & 1) return {};
                    if(deg - i / 2 <= 0) break;
                    auto ret = (*this >> i).sqrt(deg - i / 2) << (i / 2);
                    if(ret.size() < deg) ret.resize(deg, T(0));
                    return ret;
                }
            }
            return P(deg, 0);
        }

        P ret({T(1)});
        T inv2 = T(1) / T(2);
        for(int i = 1; i < deg; i <<= 1) {
            ret = (ret + pre(i << 1) * ret.inv(i << 1)) * inv2;
        }
        return ret.pre(deg);
    }

    // F(0) must be 0
    P exp(int deg = -1) const {
        assert((*this)[0] == T(0));
        const int n = (int) this->size();
        if(deg == -1) deg = n;
        P ret({T(1)});
        for(int i = 1; i < deg; i <<= 1) {
            ret = (ret * (pre(i << 1) + T(1) - ret.log(i << 1))).pre(i << 1);
        }
        return ret.pre(deg);
    }

    P pow(int64_t k, int deg = -1) const {
        const int n = (int) this->size();
        if(deg == -1) deg = n;
        for(int i = 0; i < n; i++) {
            if((*this)[i] != T(0)) {
                T rev = T(1) / (*this)[i];
                P C(*this * rev);
                P D(n - i);
                for(int j = i; j < n; j++) D[j - i] = C[j];
                D = (D.log() * k).exp() * (*this)[i].pow(k);
                P E(deg);
                if(i * k > deg) return E;
                auto S = i * k;
                for(int j = 0; j + S < deg && j < D.size(); j++) E[j + S] = D[j];
                return E;
            }
        }
        return *this;
    }


    T eval(T x) const {
        T r = 0, w = 1;
        for(auto &v : *this) {
            r += w * v;
            w *= x;
        }
        return r;
    }
};
/* how to use
// you can use ntt instead.
//with fft
ArbitraryModConvolution<mint> arb;
using FPS = FormalPowerSeries< mint >;
auto mult = [&](const FPS::P &a, const FPS::P &b) {
        auto ret = arb.multiply(a, b);
        return FPS::P(ret.begin(), ret.end());
};
FPS::set_fft(mult);
//with loop
using FPS = FormalPowerSeries< mint >;
auto mult = [&](const FPS::P &a, const FPS::P &b) {
    vi c(a.size()+b.size()-1);
    rep(i, a.size())rep(j, b.size()){
        c[i+j]+=a[i]*b[j];
    }
    return FPS::P(c.begin(), c.end())
};
FPS::set_fft(mult);
//example
FormalPowerSeries< mint > X(K + 1);
*/
//main
template< typename T >
struct Combination {
    vector< T > _fact, _rfact, _inv;

    Combination(int sz) : _fact(sz + 1), _rfact(sz + 1), _inv(sz + 1) {
        _fact[0] = _rfact[sz] = _inv[0] = 1;
        for(int i = 1; i <= sz; i++) _fact[i] = _fact[i - 1] * i;
        _rfact[sz] /= _fact[sz];
        for(int i = sz - 1; i >= 0; i--) _rfact[i] = _rfact[i + 1] * (i + 1);
        for(int i = 1; i <= sz; i++) _inv[i] = _rfact[i] * _fact[i - 1];
    }

    inline T fact(int k) const { return _fact[k]; }

    inline T rfact(int k) const { return _rfact[k]; }

    inline T inv(int k) const { return _inv[k]; }

    T P(int n, int r) const {
        if(r < 0 || n < r) return 0;
        return fact(n) * rfact(n - r);
    }

    T C(int p, int q) const {
        if(q < 0 || p < q) return 0;
        return fact(p) * rfact(q) * rfact(p - q);
    }

    T H(int n, int r) const {
        if(n < 0 || r < 0) return (0);
        return r == 0 ? 1 : C(n + r - 1, r);
    }
};
int main(){
    cin.tie(0);
    ios::sync_with_stdio(0);
    cout << setprecision(20);
    int n; cin>>n;
    Combination<mint> comb(n);
    vi p(n);
    rep(i, n){
        cin>>p[i];
        p[i]--;
    }
    UnionFind tree(n);
    rep(i, n){
        if(p[i]>=0)tree.unite(i, p[i]);
    }
    vector<bool> used(n, true);
    rep(i, n){
        if(p[i]==-2){
            used[tree.find(i)] = false;
        }
    }
    vi sizes;
    rep(i, n){
        if(tree.find(i)==i&&(not used[i])){
            sizes.pb(tree.size(i));
        }
    }
    int m = sizes.size();
    using FPS = FormalPowerSeries< mint >;
    auto mult = [&](const FPS::P &a, const FPS::P &b) {
        FPS c(a.size()+b.size()-1);
        rep(i, a.size())rep(j, b.size()){
            c[i+j]+=a[i]*b[j];
        }
        return c;
    };
    FPS::set_fft(mult);
    FPS F(1);
    F[0] = 1;
    rep(i, m){
        FPS G(2);
        G[0] = 1;
        G[1] = sizes[i];
        F*=G;
    }
    mint ans = mypow<mint>(n-1, m)*(n-1)+mypow<mint>((n-1), m);
    rep(i, 2, m+1){
        ans-=F[i]*comb.fact(i-1)*mypow<mint>((n-1), (m-i));
    }
    rep(i, m){
        ans-=(mint)(sizes[i]-1)*mypow<mint>((n-1), (m-1));
    }
    rep(i, n){
        if(tree.find(i)==i&&used[i]){
            ans-=mypow<mint>((n-1), m);
        }
    }
    cout << ans << endl;
} 

Measured memory_bytes: 3916