/+ dub.sdl:
    name "A"
    dependency "dcomp" version=">=0.4.0"
+/

import std.algorithm, std.conv, std.range, std.stdio;
// import dcomp.scanner;

// import dcomp.graph.namori;


// import dcomp.numeric.modint, dcomp.numeric.primitive, dcomp.array, dcomp.algorithm;
import std.typecons;
alias Mint = ModInt!(10^^9 + 7);
alias E = Tuple!(int, "to");
int main() {
    auto sc = new Scanner(stdin);
    int n;
    int[] a;
    sc.read(n, a); a[] -= 1;
    E[][] g = new E[][](n);
    bool[int[2]] edges;
    foreach (int i, d; a) {
        g[i] ~= E(d);
        g[d] ~= E(i);
        edges[[i, d].fixed] = true;
    }
    auto nmr = directedNamori(a);
    // non cycle, deg <= 2
    if (iota(n).filter!(i => !nmr.isCycle[i]).find!(i => g[i].length >= 3).empty == false) {
        writeln(0);
        return 0;
    }
    // cycle, deg <= 3
    if (iota(n).filter!(i => nmr.isCycle[i]).find!(i => g[i].length >= 4).empty == false) {
        writeln(0);
        return 0;
    }
    int[] co = new int[n];
    nmr.root.each!(i => co[i]++);

    int[int] naives;
    Mint ans = 1;
    foreach (cy; nmr.cycles) {
        int m = cy.length.to!int;
        auto cl = cy.map!(i => co[i]).array;
//        writeln(cy.map!(i => co[i]));
        if (cl.maximum == 1) {
            naives[m]++;
            continue;
        }
        foreach (int i, d; cl) {
            if (d == 1) continue;
            int len = 1;
            int j = (i+m-1)%m;
            while (cl[j] == 1) {
                len++;
                j = (j+m-1)%m;
            }
            //j, i
            if (d-1 > len) {
                writeln(0);
                return 0;
            }
            if (d-1 < len) {
                ans *= Mint(2);
            }
        }
    }

    auto fact = factTable!Mint(n+1);
    auto iFac = invFactTable!Mint(n+1);
    Mint C(int n, int k) {
        if (n < 0 || n < k) return Mint(0);
        return fact[n] * iFac[k] * iFac[n-k];
    }
    auto Mt = Mint(1).recurrence!((a, n) => a[n-1]*Mint(2*n-1)).take(n+1).array;
    foreach (k, co; naives) {
        //k-cycle * v 
        Mint wo = Mint(1);
        if (k >= 3 && k % 2 == 1) wo = Mint(2);
        ans *= iota(0, co+1, 2).map!(i =>
            C(co, i)*Mt[i/2]*pow(Mint(k), i/2)*pow(wo, co-i)
        ).sum;
    }
    writeln(ans);
    return 0;
}
/* IMPORT /Users/yosupo/Program/dcomp/source/dcomp/memory.d */
// module dcomp.memory;

MemoryPool _nowPool;

@property void nowPool(MemoryPool newPool) {
    _nowPool = newPool;
}
@property MemoryPool nowPool() {
    return _nowPool;
}

class MemoryPool {
    import core.memory : GC;
    import std.algorithm : max;

    static immutable Allign = 16;
    static immutable BIG = 2^^20;
    GC.BlkInfo[] blks;
    size_t e, idx, pos;
    this(size_t e) {
        this.e = e;
        idx = pos = 0;
    }
    void assign(size_t sz) {
        blks ~= GC.qalloc(sz);
    }
    void* malloc(size_t sz) {
        if (BIG < sz) {
            return GC.malloc(sz);
        }
        sz = (sz + Allign-1) / Allign * Allign;
        while (idx < blks.length && blks[idx].size < pos + sz) {
            idx++; pos = 0;
        }

        if (idx == blks.length) {
            assign(max(e, sz)); pos = 0;
        }
        pos += sz;
        assert(pos <= blks[idx].size);
        return cast(void *)(cast(byte *)(blks[idx].base) + pos - sz);
    }
    void allFree() {
        idx = pos = 0;
    }
}
/* IMPORT /Users/yosupo/Program/dcomp/source/dcomp/array.d */
// module dcomp.array;

// import dcomp.memory;

struct FastAppender(A, bool useMyPool = false) {
    @disable this(this); //this is not reference type(don't copy!)

    import std.algorithm : max;
    import std.range.primitives : ElementEncodingType;
    import core.stdc.string : memcpy;

    private alias T = ElementEncodingType!A;
    private T* _data;
    private size_t len, cap;

    this(T[]* arr) {
        _data = (*arr).ptr;
        len = (*arr).length;
        cap = (*arr).length;
    }
    void reserve(size_t nlen) {
        import core.memory : GC;
        if (nlen <= cap) return;
        
        void* nx;
        static if (useMyPool) {
            nx = nowPool.malloc(nlen * T.sizeof);
        } else {
            nx = GC.malloc(nlen * T.sizeof);            
        }

        cap = nlen;
        if (len) memcpy(nx, _data, len * T.sizeof);
        _data = cast(T*)(nx);
    }
    void opOpAssign(string op : "~")(T item) {
        if (len == cap) {
            reserve(max(4, cap*2));
        }
        _data[len++] = item;
    }
    void clear() {
        len = 0;
    }    
    T[] data() {
        return (_data) ? _data[0..len] : null;
    }
}

T[N] fixed(T, int N)(T[N] a) {return a;}
/* IMPORT /Users/yosupo/Program/dcomp/source/dcomp/graph/dfsinfo.d */
// module dcomp.graph.dfsinfo;

struct DFSInfo {
    int[] low, ord, par, vlis; //low, ord, parent, visitList
    int[][] tr; //dfs tree(directed)
    this(int n) {
        low = new int[n];
        ord = new int[n];
        par = new int[n];
        vlis = new int[n];
        tr = new int[][](n);
    }
}

DFSInfo dfsInfo(T)(T g) {
    import std.algorithm : min, each, filter;
    import std.conv : to;
    const int n = g.length.to!int;
    auto info = DFSInfo(n);
    with(info) {
        int co = 0;
        bool[] used = new bool[](n);
        void dfs(int p, int b) {
            used[p] = true;
            low[p] = ord[p] = co++;
            par[p] = b;

            bool rt = true;
            foreach (e; g[p]) {
                int d = e.to;
                if (rt && d == b) {
                    rt = false;
                    continue;
                }
                if (!used[d]) {
                    dfs(d, p);
                    low[p] = min(low[p], low[d]);
                } else {
                    low[p] = min(low[p], ord[d]);
                }
            }
        }
            
        foreach (i; 0..n) {
            if (used[i]) continue;
            dfs(i, -1);
        }
        par.filter!"a!=-1".each!((i, v) => tr[v] ~= i.to!int);
        ord.each!((i, v) => vlis[v] = i.to!int);
    }
    return info;
}
/* IMPORT /Users/yosupo/Program/dcomp/source/dcomp/numeric/primitive.d */
// module dcomp.numeric.primitive;

import std.traits;

T pow(T, U)(T x, U n) { return pow(x, n, T(1)); }
T pow(T, U)(T x, U n, T e) {
    while (n) {
        if (n & 1) e *= x;
        x *= x;
        n /= 2;
    }
    return e;
}

unittest {
    assert(pow(3, 5) == 243);
    assert(pow(3, 5, 2) == 486);
}

T lcm(T)(in T a, in T b) {
    import std.numeric : gcd;
    return a / gcd(a,b) * b;
}

unittest {
    assert(lcm(2, 4) == 4);
    assert(lcm(3, 5) == 15);
    assert(lcm(1, 1) == 1);
    assert(lcm(0, 100) == 0);
}

//a*T[0]+b*T[1]=T[2], T[2]=gcd
//todo: to binary extgcd
T[3] extGcd(T)(in T a, in T b) 
if (!isIntegral!T || isSigned!T) //unsignedはNG
{
    if (b==0) {
        return [1, 0, a];
    } else {
        auto e = extGcd(b, a%b);
        return [e[1], e[0]-a/b*e[1], e[2]];
    }
}

unittest {
    import std.numeric : gcd;
    foreach (i; 0..100) {
        foreach (j; 0..100) {
            auto e = extGcd(i, j);
            assert(e[2] == gcd(i, j));
            assert(e[0] * i + e[1] * j == e[2]);
        }
    }
}

T[] factTable(T)(size_t length) {
    import std.range : take, recurrence;
    import std.array : array;
    return T(1).recurrence!((a, n) => a[n-1]*T(n)).take(length).array;
}

// optimize
T[] invFactTable(T)(size_t length) {
    import std.range : take, recurrence;
    import std.array : array;
    return T(1).recurrence!((a, n) => a[n-1]/T(n)).take(length).array;
}

unittest {
    import std.stdio;
//     import dcomp.numeric.modint;
    alias Mint = ModInt!(10^^9 + 7);
    auto r = factTable!Mint(20);
    Mint a = 1;
    assert(r[0] == Mint(1));
    foreach (i; 1..20) {
        a *= Mint(i);
        assert(r[i] == a);
    }
    auto p = invFactTable!Mint(20);
    foreach (i; 1..20) {
        assert((r[i]*p[i]).v == 1);
    }
}
/* IMPORT /Users/yosupo/Program/dcomp/source/dcomp/scanner.d */
// module dcomp.scanner;

class Scanner {
    import std.stdio : File;
    import std.conv : to;
    import std.range : front, popFront, array, ElementType;
    import std.array : split;
    import std.traits : isSomeChar, isStaticArray, isArray; 
    import std.algorithm : map;
    File f;
    this(File f) {
        this.f = f;
    }
    string[] buf;
    private bool succ() {
        while (!buf.length) {
            if (f.eof) return false;
            buf = f.readln.split;
        }
        return true;
    }
    private bool readSingle(T)(ref T x) {
        if (!succ()) return false;
        static if (isArray!T) {
            alias E = ElementType!T;
            static if (isSomeChar!E) {
                //string or char[10] etc
                x = buf.front;
                buf.popFront;
            } else {
                static if (isStaticArray!T) {
                    //static
                    assert(buf.length == T.length);
                }
                x = buf.map!(to!E).array;
                buf.length = 0;                
            }
        } else {
            x = buf.front.to!T;
            buf.popFront;            
        }
        return true;
    }
    int read(T, Args...)(ref T x, auto ref Args args) {
        if (!readSingle(x)) return 0;
        static if (args.length == 0) {
            return 1;
        } else {
            return 1 + read(args);
        }
    }
}

unittest {
    import std.path : buildPath;
    import std.file : tempDir;
    import std.algorithm : equal;
    import std.stdio : File;
    string fileName = buildPath(tempDir, "kyuridenanmaida.txt");
    auto fout = File(fileName, "w");
    fout.writeln("1 2 3");
    fout.writeln("ab cde");
    fout.writeln("1.0 1.0 2.0");
    fout.close;
    Scanner sc = new Scanner(File(fileName, "r"));
    int a;
    int[2] b;
    char[2] c;
    string d;
    double e;
    double[] f;
    sc.read(a, b, c, d, e, f);
    assert(a == 1);
    assert(equal(b[], [2, 3]));
    assert(equal(c[], "ab"));
    assert(equal(d, "cde"));
    assert(e == 1.0);
    assert(equal(f, [1.0, 2.0]));
}
/* IMPORT /Users/yosupo/Program/dcomp/source/dcomp/algorithm.d */
// module dcomp.algorithm;

//[0,0,0,...,1,1,1]で、初めて1となる場所を探す。pred(l) == 0, pred(r) == 1と仮定
T binSearch(alias pred, T)(T l, T r) {
    while (r-l > 1) {
        T md = (l+r)/2;
        if (!pred(md)) l = md;
        else r = md;
    }
    return r;
}

import std.range.primitives;

Rotator!Range rotator(Range)(Range r)
if (isForwardRange!Range && hasLength!Range) {
    return typeof(return)(r);
}

struct Rotator(Range)
if (isForwardRange!Range && hasLength!Range) {
    size_t cnt;
    Range start, now;
    this(Range r) {
        cnt = 0;
        start = r.save;
        now = r.save;
    }
    this(this) {
        start = start.save;
        now = now.save;
    }
    @property bool empty() {
        return now.empty;
    }
    @property auto front() {
        assert(!now.empty);
        import std.range : take, chain;
        return chain(now, start.take(cnt));
    }
    @property Rotator!Range save() {
        return this;
    }
    void popFront() {
        cnt++;
        now.popFront;
    }
}


E minimum(alias pred = "a < b", Range, E = ElementType!Range)(Range range, E seed)
if (isInputRange!Range && !isInfinite!Range && !is(CommonType!(ElementType!Range, E) == void)) {
    import std.functional;
    while (!range.empty) {
        if (binaryFun!pred(range.front, seed)) {
            seed = range.front;
        }
        range.popFront;
    }
    return seed;
}

//should be use reduce?
ElementType!Range minimum(alias pred = "a < b", Range)(Range range)
if (isInputRange!Range && !isInfinite!Range) {
    assert(!range.empty, "range must not empty");
    auto e = range.front; range.popFront;
    return minimum!pred(range, e);
}

E maximum(alias pred = "a < b", Range, E = ElementType!Range)(Range range, E seed)
if (isInputRange!Range && !isInfinite!Range && !is(CommonType!(ElementType!Range, E) == void)) {
    import std.functional;
    while (!range.empty) {
        if (binaryFun!pred(seed, range.front)) {
            seed = range.front;
        }
        range.popFront;
    }
    return seed;
}

ElementType!Range maximum(alias pred = "a < b", Range)(Range range)
if (isInputRange!Range && !isInfinite!Range) {
    assert(!range.empty, "range must not empty");
    auto e = range.front; range.popFront;
    return maximum!pred(range, e);
}

unittest {
    assert(minimum([2, 1, 3]) == 1);
    assert(minimum!"a > b"([2, 1, 3]) == 3);
    assert(minimum([2, 1, 3], -1) == -1);
    assert(minimum!"a > b"([2, 1, 3], 100) == 100);

    assert(maximum([2, 1, 3]) == 3);
    assert(maximum!"a > b"([2, 1, 3]) == 1);
    assert(maximum([2, 1, 3], 100) == 100);
    assert(maximum!"a > b"([2, 1, 3], -1) == -1);
}

bool[ElementType!Range] toMap(Range)(Range r) {
    import std.algorithm : each;
    bool[ElementType!Range] res;
    r.each!(a => res[a] = true);
    return res;
}
/* IMPORT /Users/yosupo/Program/dcomp/source/dcomp/graph/namori.d */
// module dcomp.graph.namori;

// import dcomp.graph.dfsinfo;

struct Namori {
    bool[] isCycle;
    int[][] cycles;
    int[] root;
    this(int n) {
        isCycle = new bool[n];
        root = new int[n];
    }
}

Namori namori(T)(T g) {
    import std.algorithm : find, each;
    import std.range;
    import std.conv : to;
    int n = g.length.to!int;
    auto info = dfsInfo(g);
    auto nmr = Namori(n);
    with (nmr) {
        //find self loop
        foreach (i; 0..n) {
            if (g[i].find!(e => e.to == i).empty) continue;
            isCycle[i] = true;
            cycles ~= [i];
        }
        foreach (p; info.vlis) {
            if (info.low[p] == info.ord[p]) continue;
            if (g[p].length == info.tr[p].length+1) continue;
            int[] v;
            int nw = p;
            while (info.ord[nw] != info.low[p]) {
                v ~= nw;
                nw = info.par[nw];
            }
            v ~= nw;
            v.each!(i => isCycle[i] = true);
            cycles ~= v;
        }
        bool[] used = new bool[n];
        void dfs(int p, int b, int r) {
            if (used[p]) return;
            used[p] = true;
            root[p] = r;
            foreach (e; g[p]) {
                int d = e.to;
                if (d == b) continue;
                if (isCycle[d]) continue;
                dfs(d, p, r);
            }
        }
        foreach (i; 0..n) {
            if (!isCycle[i]) continue;
            dfs(i, -1, i);
        }
    }
    return nmr;
}

Namori directedNamori(int[] g) {
    import std.algorithm : find, each;
    import std.range;
    import std.conv : to;
    int n = g.length.to!int;
    auto nmr = Namori(n);
    with (nmr) {
        int[] used = new int[n]; used[] = -1;
        foreach (i; 0..n) {
            if (used[i] != -1) continue;
            int j = i;
            while (used[j] == -1) {
                used[j] = i;
                j = g[j];
            }
            if (used[j] != i) continue;
            int k = j;
            int[] cy = [];
            while (true) {
                cy ~= k;
                isCycle[k] = true;
                if (g[k] == j) break;
                k = g[k];
            }
            cycles ~= cy;
        }
        int[] dp = new int[n]; dp[] = -1;
        int dfs(int p) {
            if (isCycle[p]) return p;
            if (dp[p] != -1) return dp[p];
            return dp[p] = dfs(g[p]);
        }
        foreach (i; 0..n) {
            root[i] = dfs(i);
        }
    }
    return nmr;
}
/* IMPORT /Users/yosupo/Program/dcomp/source/dcomp/numeric/modint.d */
// module dcomp.numeric.modint;

// import dcomp.numeric.primitive;

struct ModInt(uint MD) if (MD < int.max) {
    import std.conv : to;
    uint v;
    this(int v) {this.v = (long(v)%MD+MD)%MD;}
    this(long v) {this.v = (v%MD+MD)%MD;}
    auto normS(uint x) {return (x<MD)?x:x-MD;}
    auto make(uint x) {ModInt m; m.v = x; return m;}
    auto opBinary(string op:"+")(ModInt r) {return make(normS(v+r.v));}
    auto opBinary(string op:"-")(ModInt r) {return make(normS(v+MD-r.v));}
    auto opBinary(string op:"*")(ModInt r) {return make( (long(v)*r.v%MD).to!uint );}
    auto opBinary(string op:"/")(ModInt r) {return this*inv(r);}
    auto opOpAssign(string op)(ModInt r) {return mixin ("this=this"~op~"r");}
    static ModInt inv(ModInt x) {return ModInt(extGcd!int(x.v, MD)[0]);}
    string toString() {return v.to!string;}
}

unittest {
    static assert( is(ModInt!(uint(1000000000) * 2))); //not overflow
    static assert(!is(ModInt!(uint(1145141919) * 2))); //overflow!
    alias Mint = ModInt!(10^^9+7);
    // negative check
    assert(Mint(-1).v == 10^^9 + 6);
    assert(Mint(-1L).v == 10^^9 + 6);

    Mint a = 48;
    Mint b = Mint.inv(a);
    assert(b.v == 520833337);

    Mint c = Mint(15);
    Mint d = Mint(3);
    assert((c/d).v == 5);
}

struct DModInt {
    import std.conv : to;
    uint MD, v;
    this(int v, uint md) {
        MD = md;
        this.v = ((long(v)%MD+MD)%MD).to!uint;
    }
    this(long v, uint md) {
        MD = md;
        this.v = ((v%MD+MD)%MD).to!uint;
    }
    auto normS(uint x) {return (x<MD)?x:x-MD;}
    auto make(uint x) {DModInt m; m.MD = MD; m.v = x; return m;}
    auto opBinary(string op:"+")(DModInt r) {
        assert(MD == r.MD);
        return make(normS(v+r.v));
    }
    auto opBinary(string op:"-")(DModInt r) {
        assert(MD == r.MD);
        return make(normS(v+MD-r.v));
    }
    auto opBinary(string op:"*")(DModInt r) {
        assert(MD == r.MD);
        return make((long(v)*r.v%MD).to!uint);
    }
    auto opBinary(string op:"/")(DModInt r) {
        assert(MD == r.MD);
        return this*inv(r);
    }
    auto opOpAssign(string op)(DModInt r) {return mixin ("this=this"~op~"r");}
    static DModInt inv(DModInt x) {
        return DModInt(extGcd!int(x.v, x.MD)[0], x.MD);
    }
    string toString() {return v.to!string;}
}

unittest {
    immutable MD = 10^^9 + 7;
    alias Mint = DModInt;
    //negative check
    assert(Mint(-1, MD).v == 10^^9 + 6);
    assert(Mint(-1L, MD).v == 10^^9 + 6);
    Mint a = Mint(48, MD);
    Mint b = Mint.inv(a);
    assert(b.v == 520833337);
    Mint c = Mint(15, MD);
    Mint d = Mint(3, MD);
    assert((c/d).v == 5);
}


Measured memory_bytes: 16060