We use many equations below. If you want an intuitive understanding, see Regarding this problem as a graph problem.

By symmetry, one can find the sum of \(S_1\) over \(N^N\) candidates of \(A\) and multiply it by \(N\).

Let \(L\) be the maximum \(j\) such that \(B_{1,1},B_{1,2},\dots\) and \(B_{1,j}\) are pairwise distinct. We exhaustively inspect \(L\) between \(1\) and \(N\), inclusive.

When \(L=l\), there are \({}_{N-1}P_{l-1}\) ways to choose \(B_{1,2},\dots,B_{1,l}\). Moreover, \(B_{1,l+1}\) coincides with one of \(B_{1,1},B_{1,2},\dots,B_{1,l}\). If \(B_{1,l+1}=B_{1,p}\ (1\leq p \leq l)\), by definition of \(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},\]

so \(B_{1,p},\dots\) and \(B_{1,l}\) occurs repeatedly as:

\[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).\]

Therefore, \(B_{1,1},B_{1,2},\dots\), and \(B_{1,p-1}\) are the integers that occurs exactly once in \(B_1\), and \(S_1=p-1\).

Now let us fix \(B_{1,1},B_{1,2},\dots,B_{1,l+1}\) and count the number of corresponding \(A\). By the discussion above, \(A_{B_{1,1}},A_{B_{1,2}},\dots,A_{B_{1,l}}\) is already determined. On the other hand, any of the other \((N-l)\) elements do not affect \(B_1\) and \(S_1\).

Therefore, the sum of \(S_1\) over all \(A\) such that \(A=l\) is

\[{}_{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}.\]

Since \({}_{N-1}P_{l-1}\) over \(l=1,2,...,N\) can be enumerated in a total of \(O(N)\) time, the problem can be solved in a total of \(O(N)\) time.

[Bonus] Regarding this problem as a graph problem

Consider a graph \(G\) with \(N\) vertices and \(N\) edges, whose \(i\)-th edge goes from vertices \(i\) to \(A_i\). As shown in the graph below, each connected component of \(G\) consists of a cycle and some rooted trees rooted at vertices on the cycle (such a graph is called “functional graph”, or informally “Namori graph” in Japan).

Let us interpret our problem on this graph \(G\). Then \(B_i\) is “the sequence of the indices of the vertices to visit when you start from vertex \(i\) and repeat traversing the outgoing edge \((10^{100}-1)\) times.” During the travel, once you reached the cycle, you repeatedly go around the cycle, so \(S_i\) appears to be “the number of vertices you visit before you reach the cycle for the first time if you start from vertex \(i\) and repeatedly traverse the outgoing edge.” Since \(B_i\) only contains the vertices on “the path” (from vertex \(i\) to the cycle) and the cycle, all that left is an appropriate combinatoric computation (specifically, choosing the vertices on the path and the cycle, and then determining the destination of vertices not contained in the path nor the cycle) to find the sum of \(S_i\). This approach ends up in the equations described in the main text in this editorial. (\(L\) in the main text corresponds to “the sum of the lengths of the path and the cycle” in this part, and \(p\) “the length of the path.”

Sample code (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;
}