この問題はダブリングにより解くことができます。

端末 \(S\) から \(K\) 回操作を行ったときに得られる値を \(f(K,S)\) と表すことにします。また、端末 \(S\) から \(K\) 回移動した後の端末の番号を \(g(K,S)\) とします。

操作を \(n\) 回目までと \(n\) 回目よりあとに分けて考えることで、\(n\leq K\) を満たす任意の \(n\) について、\(f(K,S)=f(n,S)\times 10^{K-n}+f(K-n, g(n,S))\) が成り立つことがわかります。

（\(n\) 回目の操作が終わった時点で端末 \(g(n,S)\) にいるため、そこから残りの \(K-n\) 回の操作を行って得られるのは \(f(K-n,g(n,S))\) です。これが \(n\) 回目までの操作で得られている \(f(n,S)\) の後ろに \(K-n\) 桁つくので、全体では先程の式になります）

予め全ての \(S\) と全ての \(0 \leq i \leq B\) に対して \(f(2^i,S)\) と \(g(2^i,S)\) を前計算します。

この前計算は

• \(g(2^i,S)=g(2^{i-1},g(2^{i-1},S))\)
• \(f(2^i,S)=f(2^{i-1},S)\times 10^{2^{i-1}}+f(2^{i-1},g(2^{i-1},S))\)

の漸化式により行えるため、(\(10^{2^{i-1}}\) も適切に計算すると) 全体で \(O(BN)\) で行えます。

この前計算の下、 \(K\) を \(2\) ベキに分解することで \(f(K,S)\) を \(O(\log K)\) 時間で計算することができます。

\(B=\log_2 \max K_i\) と取ることで、全体で \(O((N+Q)\log\max K_i)\) 時間で解くことができます。

実装例 (C++)

#include<bits/stdc++.h>
#include<atcoder/modint>
using namespace std;
using mint = atcoder::modint;

int main(){
  int n, q, m;
  cin >> n >> q >> m;
  mint::set_mod(m);
  vector<int>d(n), p(n);
  for(int i=0; i<n; i++){
    cin >> d[i] >> p[i];
    p[i]--;
  }

  vector<vector<mint>>f(61, vector<mint>(n));
  vector<vector<int>>g(61, vector<int>(n));
  vector<mint>pow10(61);  // pow10[i] = 10^2^i
  for(int i=0; i<n; i++){
    f[0][i] = d[i];
    g[0][i] = p[i];
  }
  pow10[0] = 10;
  
  for(int k=1; k<=60; k++){
    for(int i=0; i<n; i++){
      f[k][i] = f[k-1][i] * pow10[k-1] + f[k-1][g[k-1][i]];
      g[k][i] = g[k-1][g[k-1][i]];
    }
    pow10[k] = pow10[k-1] * pow10[k-1];
  }
  
  for(int i=0; i<q; i++){
    int s, k;
    cin >> s >> k;
    s--;
    mint ans = 0;
    for(int kk=0; kk<=60; kk++){
      if(k & (1<<kk)){
        ans = ans * pow10[k] + f[k][s];
        s = g[k][s];
      }
    }
    cout << ans.val() << endl;
  }
}

実装例 (Python)

N, Q, M = map(int, input().split())
D, P = [], []
for _ in range(N):
  d, p = map(int, input().split())
  D.append(d)
  P.append(p-1)

f = [[0]*N for _ in range(61)]
g = [[0]*N for _ in range(61)]
pow10 = [0] * 61  # pow10[i] = pow(10, 2**i, M)

f[0] = D
g[0] = P
pow10[0] = 10

for k in range(1, 61):
  for i in range(N):
    f[k][i] = (f[k-1][i] * pow10[k-1] + f[k-1][g[k-1][i]]) % M
    g[k][i] = g[k-1][g[k-1][i]]
  pow10[k] = pow10[k-1]**2 % M

for _ in range(Q):
  S, K = map(int, input().split())
  S -= 1
  ans = 0
  for k in range(61):
    if K & (1<<k):
      ans = (ans * pow10[k] + f[k][S]) % M
      S = g[k][S]
  print(ans)