AOJ 2606 : Perm Query

順列はいくつかのサイクルで表記できる。
[l, r)に含まれるサイクルの長さをc[l],...,c[r-1]のように表すと
lcm(c[l], ..., c[r])回だけ、各要素は値が変わる。i番目の要素は
lcm/c[i] 回だけサイクルを回るので、そのサイクルに含まれる数の和をS[i]とすると
i番目の要素についての和は
S[i]*lcm/c[i]
l≦i＜rで和を取って
Σ[l≦i＜r](S[i]*lcm/c[i])
がクエリ[l, r)に対する解。
問題はlcmの求め方で、lcmは大きい数になりうるのでmod 10^9+7で計算したい。
つまり
[l, r)に含まれるすべてのサイクルの長さについてlcmを計算したい。
[l, r)はO(N)の長さがあって1要素ずつ見ていくと間に合いそうにない。
そこでサイクルの長さごとに見ていく。
サイクルの長さの種類はたかだかO(√N)個しかない。なぜなら、サイクルの長さの和はNなので、
サイクルの長さの種類が最大となるのは
1+2+...+(k-1)+k = N
のような場合であるため。
よって、O(√N)個のすべてのサイクルの長さについて、そのサイクルの出現数を部分和で持っておけば、
各サイクルの長さが[l, r)で出現するかどうかすぐわかる。
事前にサイクルの長さに対して素因数分解しておけば、素因数の数はたかだかO(log(N))個なので
各クエリに対してO(log(N)√N)
全体でO(Qlog(N)√N)

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    int N, Q;
    cin >> N >> Q;
    vi p(N), vis(N);
    // S[i]/c[i]の部分和
    vector<mint> sm(N + 1);
    // 素因数分解の結果
    unordered_map<int, map<int, int> > fa;
    // lns[i]は長さiのサイクルの出現の部分和
    unordered_map<int, vi> lns;
    rep(i, N)cin >> p[i], p[i]--;
    rep(i, N)if (!vis[i]) { // 新しいサイクルを見つけた
        vi a;
        ll s = 0;
        int x = i;
        do {
            a.push_back(x);
            vis[x] = 1;
            // サイクルに含まれる数の和
            s += x + 1;
            x = p[x];
        } while (x != i);
        int len = sz(a);
        if (!fa.count(len)) {
            fa[len] = factorize(len);
            lns[len] = vi(N + 1);
        }
        each(y, a) {
            // S[i]/c[i]
            sm[y + 1] = mint(s) / len;
            // 長さlenのサイクルはy番目を含む
            lns[len][y + 1] = 1;
        }
    }
    // 部分和の計算
    rep(i, N)sm[i + 1] += sm[i];
    each(q, lns) {
        vi &v = q.second;
        // 出現数の部分和
        rep(i, N)v[i + 1] += v[i];
    }
    while (Q--) {
        int l, r;
        cin >> l >> r;
        --l;
        mint s = sm[r] - sm[l], lcm = 1;
        // 各素因数について、素因数分解したときに最も多く現れたときの指数
        unordered_map<int, int> ma;
        each(q, lns) {
            int len = q.first;
            vi &v = q.second;
            if (v[r] - v[l]) {
                each(ii, fa[len]) {
                    // ここが最も時間がかかる
                    smax(ma[ii.first], ii.second);
                }
            }
        }
        // 素因数分解→lcm
        each(q, ma) while (q.second--)lcm *= q.first;
        // s=Σ(S[i]/c[i])
        cout << s*lcm << endl;
    }
}