以下、数式を多く用いて解説します。先に感覚的な理解を得たい方は、本問題をグラフの問題として捉える方法 をご参照ください。

対称性より、\(N^N\) 通りの \(A\) に対する \(S_1\) の総和を求めた上で、それを \(N\) 倍すれば良いです。

「\(B_{1,1},B_{1,2},\dots,B_{1,j}\) に含まれる数が全て相異なるような最大の \(j\)」を \(L\) と置きます。\(L\) を \(1\) 以上 \(N\) 以下の範囲で全探索します。

\(L=l\) のとき、\(B_{1,2},\dots,B_{1,l}\) として有り得るものは \({}_{N-1}P_{l-1}\) 通りあります。また、\(B_{1,l+1}\) は \(B_{1,1},B_{1,2},\dots,B_{1,l}\) のいずれかと一致しますが、\(B_{1,l+1}=B_{1,p}\ (1\leq p \leq l)\) のとき、\(B_1\) の定義より

\[A_{B_{1,1}}=B_{1,2},A_{B_{1,2}}=B_{1,3},\dots,A_{B_{1,l-1}}=B_{1,l}\]

\[A_{B_{1,l}}=B_{1,p}\]

ですから、

\[B_1=(B_{1,1},\dots,B_{1,p-1},B_{1,p},\dots,B_{1,l},B_{1,p},\dots,B_{1,l},B_{1,p},\dots,B_{1,l},\dots)\]

というように、\(B_{1,p},\dots,B_{1,l}\) のかたまりが繰り返されることになります。よって、\(B_1\) のなかでちょうど \(1\) 度だけ出てくる数は \(B_{1,1},B_{1,2},\dots,B_{1,p-1}\) であり、\(S_1=p-1\) です。

では、\(B_{1,1},B_{1,2},\dots,B_{1,l+1}\) を \(1\) つ固定したとき、対応する \(A\) が何通りあるか考えてみましょう。上記の議論より、\(A_{B_{1,1}},A_{B_{1,2}},\dots,A_{B_{1,l}}\) については既に確定しています。逆に、それ以外の \(N-l\) 個の要素が何であっても\(B_1\) および \(S_1\) には影響しません。

よって、\(L=l\) であるような \(A\) に対する \(S_1\) の総和は、

\[{}_{N-1}P_{l-1}N^{N-l}\sum_{p=1}^{l}(p-1)={}_{N-1}P_{l-1}N^{N-l}\times\frac{l(l-1)}{2}\]

です。\(l=1,2,...,N\) に対する \({}_{N-1}P_{l-1}\) は \(O(N)\) で列挙できるため、本問題を \(O(N)\) で解くことができました。

<参考> 本問題をグラフの問題として捉える方法

\(N\) 頂点 \(N\) 辺のグラフであって、\(i\) 番目の辺が頂点 \(i\) から \(A_i\) に伸びているようなグラフ\(G\) を考えます。\(G\) は以下の図のように、各連結成分が \(1\) つのサイクルおよびサイクル上の頂点を根としたいくつかの根付き木からなるようなグラフになります (このようなグラフは Functional Graph と呼ばれ、俗になもりグラフとも呼ばれています)。

本問題の設定をこのグラフ \(G\) 上で解釈してみましょう。すると、\(B_i\) は「頂点 \(i\) から始めて、今いる頂点から出ている辺の方向に進むことを \(10^{100}-1\) 回繰り返したとき、通過した頂点の番号を通過した順番に並べた数列」であると言えます。このとき、一度サイクルに到達した後はサイクル上を繰り返し回り続けますから、\(S_i\) は「頂点 \(i\) から始めて辺を辿って行ったとき、最初にサイクルに到達するまでに通る頂点の数」であることが分かります。\(B_i\) に含まれる頂点は「頂点 \(i\) からサイクルまでのパス (以降単にパスと呼びます)」およびサイクル上の頂点のみですから、パスとサイクルの長ささえ固定してしまえば、あとは適切な組合せ計算 (具体的には、パスとサイクル上にある頂点を選ぶ + パスにもサイクルにも含まれない頂点から伸びる辺の行き先を決める) によって \(S_i\) の総和が計算できそうだという見通しが立ちます。この方針を数式で整理していくと、本解説の本文のようになります (本文中の \(L\) は補足中の「パスとサイクルの長さの合計」に、本文中の \(p\) は補足中の「パスの長さ」にそれぞれ対応しています)。

実装例 (C++) :

#include<bits/stdc++.h>
#include<atcoder/modint>

using namespace std;
using namespace atcoder;

using ll = long long;
using mint = modint;

int main() {
    int n, m;
    cin >> n >> m;
    mint::set_mod(m);
    
    vector<mint> pow_n(n, 1);
    for (int i = 0; i < n - 1; i++) {
        pow_n[i + 1] = pow_n[i] * n;
    }
    
    mint ans = 0;
    mint perm = 1;
    for (int l = 1; l <= n; l++) {
        ans += perm * pow_n[n - l] * ((ll) l * (l - 1) / 2);
        perm *= n - l;
    }
    ans *= n;
    cout << ans.val() << endl;
}