#include<iostream>
#include<cstring>
#include<vector>
#include<array>
#include<set>
#include<algorithm>
using namespace std;
using LL = long long;
const int mod = 1e9 + 7;

void add(int &a, int b){
    a += b;
    if (a >= mod) a -= mod;
}

int mul(int a, int b){
    return 1LL * a * b % mod;
}

const int maxn = 5e5 + 5;
vector<int> factor[maxn];
int inv[maxn];

int main(){

#ifdef LOCAL
    freopen("data.in", "r", stdin);
    freopen("data.out", "w", stdout);
#endif

    cin.tie(0);
    cout.tie(0);
    ios::sync_with_stdio(0);

    inv[1] = 1;
    for(int i = 2; i < maxn; i++){
        inv[i] = mul(mod - mod / i, inv[mod % i]);
    }

    for(int i = 2; i < maxn; i++){
        if (factor[i].empty()){
            for(int j = i; j < maxn; j += i){
                factor[j].push_back(i);
            }
        }
    }

    int T;
    cin >> T;
    while(T--){
        int n;
        cin >> n;
        vector<int> p(n + 1);
        for(int i = 1; i <= n; i++) cin >> p[i];
        vector<bool> v(n + 1);
        vector<int> cnt(n + 1);
        vector<array<int, 3> > mx(n + 1);
        for(int i = 1; i <= n; i++){
            if (!v[i]){
                int len = 0;
                for(int j = i; !v[j]; j = p[j]){
                    len += 1;
                    v[j] = true;
                }
                cnt[len] += 1;
                int t = len;
                for(auto x : factor[len]){
                    int c = 0;
                    while(t % x == 0) t /= x, c += 1;
                    for(int j = 0; j < 3; j++){
                        if (c > mx[x][j]){
                            swap(c, mx[x][j]);
                        }
                    }
                }
            }
        }
        vector<int> cand;
        for(int i = 1; i <= n; i++){
            if (cnt[i]){
                cand.push_back(i);
            }
        }
        int lcm = 1;
        for(int i = 1; i <= n; i++){
            for(int j = 0; j < mx[i][0]; j++){
                lcm = mul(lcm, i);
            }
        }
        int ans = 0;
        for(int i = 0; i < cand.size(); i++){
            for(int j = i; j < cand.size(); j++){
                int x = cand[i], y = cand[j];
                int tot = 0;
                if (x == y) tot = mul(mul(cnt[x], cnt[x] - 1), mul(x, x));
                else tot = mul(2, mul(mul(cnt[x], x), mul(cnt[y], y)));

                if (tot == 0) continue;

                auto p = factor[x];
                p.insert(p.end(), factor[y].begin(), factor[y].end());
                p.insert(p.end(), factor[x + y].begin(), factor[x + y].end());
                sort(p.begin(), p.end());
                p.erase(unique(p.begin(), p.end()), p.end());
                
                int z = x + y;

                int nlcm = lcm;

                for(auto u : p){
                    multiset<int> s(mx[u].begin(), mx[u].end());
                    multiset<int>::iterator it;
                    for(auto v : {x, y}){
                        int c = 0;
                        while(v % u == 0) v /= u, c += 1;
                        it = s.find(c);
                        if (it != s.end()) s.erase(it);
                    }
                    {
                        int c = 0;
                        while(z % u == 0) z /= u, c += 1;
                        s.insert(c);
                    }
                    int cur = *s.rbegin();
                    if (cur < mx[u][0]){
                        for(int i = 0; i < mx[u][0] - cur; i++){
                            nlcm = mul(nlcm, inv[u]);
                        }
                    }
                    else{
                        for(int i = 0; i < cur - mx[u][0]; i++){
                            nlcm = mul(nlcm, u);
                        }
                    }
                }
                add(ans, mul(nlcm, tot));
            }
        }
        cout << ans << '\n';
    }

}