#include <iostream>
#include <vector>
#include <algorithm>
#include <string>
#include <set>

#define REP(i, a, b) for (int i = int(a); i < int(b); i++)
#ifdef _DEBUG_
#define dump(val) cerr << __LINE__ << ":\t" << #val << " = " << (val) << endl
#else
#define dump(val)
#endif

using namespace std;

typedef long long int ll;

template<typename T>
vector<T> make_v(size_t a, T b) {
    return vector<T>(a, b);
}
template<typename... Ts>
auto make_v(size_t a, Ts... ts) {
    return vector<decltype(make_v(ts...))>(a, make_v(ts...));
}

class UnionFind {
public:
    vector<int> par;
    int N;
    UnionFind(int n) {
        par.resize(n, -1);
        N = n;
    }
    int Find(int n) {
        return par[n] < 0 ? n : par[n] = Find(par[n]);
    }
    bool Union(int x, int y) {
        x = Find(x);
        y = Find(y);
        if (x == y) return false;
        if (-par[x] < -par[y]) swap(x, y);
        par[y] += par[x];
        par[x] = y;
        N--;
        return true;
    }
    bool Same(int x, int y) {
        return Find(x) == Find(y);
    }
    int size(int x) {
        return -par[Find(x)];
    }
    int size() {
        return N;
    }
};

template<ll MOD = 1000000007>
class ModInt {
    ll n;
    static ModInt pow(ModInt x, ll p) {
        if (p == 0) {
            return 1;
        } else if (p % 2) {
            return x * pow(x, p - 1);
        } else {
            auto t = ModInt::pow(x, p / 2);
            return t * t;
        }
    }
    ModInt reverse() const {
        return ModInt::pow(*this, MOD - 2);
    }

public:
    ModInt()
        : n(0) {}
    ModInt(ll _n)
        : n(_n % MOD) {}
    ModInt operator+=(const ModInt &m) {
        n += m.n;
        if (n >= MOD) {
            n -= MOD;
        }
        return *this;
    }
    ModInt operator-=(const ModInt &m) {
        n -= m.n;
        if (n < 0) {
            n += MOD;
        }
        return *this;
    }
    ModInt operator*=(const ModInt &m) {
        n *= m.n;
        if (n >= MOD) {
            n %= MOD;
        }
        return *this;
    }
    ModInt operator/=(const ModInt &m) {
        return n = (m.reverse() * n).n;
    }

    ModInt operator+=(const ll l) {
        return (*this) += ModInt(l);
    }
    ModInt operator-=(const ll l) {
        return (*this) -= ModInt(l);
    }
    ModInt operator*=(const ll l) {
        return (*this) *= ModInt(l);
    }
    ModInt operator/=(const ll l) {
        return (*this) /= ModInt(l);
    }
    ModInt operator+(const ModInt &m) {
        auto t = *this;
        return t += m;
    }
    ModInt operator-(const ModInt &m) {
        auto t = *this;
        return t -= m;
    }
    ModInt operator*(const ModInt &m) {
        auto t = *this;
        return t *= m;
    }
    ModInt operator/(const ModInt &m) {
        auto t = *this;
        return t /= m;
    }
    ModInt operator+(const ll l) {
        auto t = *this;
        return t += l;
    }
    ModInt operator-(const ll l) {
        auto t = *this;
        return t -= l;
    }
    ModInt operator*(const ll l) {
        auto t = *this;
        return t *= l;
    }
    ModInt operator/(const ll l) {
        auto t = *this;
        return t /= l;
    }
    ModInt operator=(const ll l) {
        n = l % MOD;
        if (n < 0) {
            n += MOD;
        }
        return *this;
    }
    friend ostream &operator<<(ostream &out, const ModInt &m) {
        out << m.n;
        return out;
    }
    friend istream &operator>>(istream &in, ModInt &m) {
        ll l;
        in >> l;
        m = l;
        return in;
    }
};

using mint = ModInt<>;

mint fast_pow(mint x, int n) {
    if (n == 0) {
        return 1;
    } else if (n % 2) {
        return fast_pow(x, n - 1) * x;
    } else {
        x = fast_pow(x, n / 2);
        return x * x;
    }
}

int main() {
    cin.tie(nullptr);
    ios::sync_with_stdio(false);
    int N;
    while (cin >> N, N) {
        vector<int> A(N);
        REP(i, 0, N) {
            cin >> A[i];
            A[i]--;
        }
        bool ok = true;
        {
            vector<set<int>> V(N);
            REP(i, 0, N) {
                V[i].insert(A[i]);
                V[A[i]].insert(i);
            }
            REP(i, 0, N) {
                ok &= (V[i].size() <= 2);
            }
        }

        UnionFind uf(N);
        REP(i, 0, N) {
            if (A[A[i]] == i || !uf.Same(A[i], i)) {
                uf.Union(i, A[i]);
            } else {
                ok = false;
            }
        }
        if (!ok) {
            cout << 0 << endl;
            continue;
        }

        int p = 0, q = 0;
        {
            vector<bool> V(N, false);
            REP(i, 0, N) {
                int par = uf.Find(i);
                if (!V[par]) {
                    V[par] = true;
                    if (uf.size(par) == 2) {
                        p++;
                    } else {
                        q++;
                    }
                }
            }
        }

        vector<mint> factor(N + 1, 1);
        REP(i, 0, N) {
            factor[i + 1] *= factor[i] * (i + 1);
        }

        auto nCr = [&](int n, int r) {
            return factor[n] / (factor[r] * factor[n - r]);
        };

        mint ans = 0;
        if (p == 0) {
            ans = fast_pow(2, q) * factor[q + 1];
        } else {
            if (q == 0) {
                ans = fast_pow(2, p) * factor[p] * (p / 2 + 1);
            } else {
                if (q >= 2) {
                    REP(i, 0, p + 1) {
                        ans += nCr(p + q - i - 2, p - i) * fast_pow(2, p + q) * factor[p] * factor[q] * (i / 2 + 1) * (p + q - i - 1);
                    }
                }
                if (q >= 1) {
                    REP(i, 0, p + 1) {
                        ans += nCr(p + q - i - 1, p - i) * fast_pow(2, p + q) * factor[p] * factor[q] * (i / 2 + 1) * 2;
                    }
                }
            }
        }
        cout << ans << endl;
    }
    return 0;
}



Measured memory_bytes: 15256