Let \(c ^ k _ {i,j}\) denote the total cost to travel from vertex \(i\) to vertex \(j\) via exactly \(k\) edges. Consider the sequence \(c ^ k=(c ^ k _ {i,j}) _ {1\le i\le N,1\le j\le N}\) indexed by \((i,j)\). We have \(c ^ 1=(C _ {i,j}) _ {1\le i\le N,1\le j\le N}\).

Given \(c ^ a\) and \(c ^ b\), the algorithm yields \(c ^ {a+b}\):

• Initialize \(x _ {i,j}\leftarrow\infty\).
• For all \((i,j,k)\in\lbrace1,2,\ldots,N\rbrace ^ 3\), repeat the following:
  • Set \(x _ {i,k}\leftarrow\min\lbrace x _ {i,k},c ^ a _ {i,j}+c ^ b _ {j,k}\rbrace\).
• \(x\) is the sought \(c ^ {a+b}\).

This can be interpreted as a matrix product on the min-plus semiring.

This algorithm enables to find the answer using fast exponentiation in \(O(N ^ 3\log K)\) time.

The following is sample code:

#include <iostream>
#include <vector>
#include <ranges>
#include <algorithm>

int main() {
    using namespace std;
    unsigned N, K;
    cin >> N >> K;
    vector C(N, vector<unsigned long>(N));

    for (auto&& row : C)
        for (auto&& c : row)
            cin >> c;

    // Finds results of (a + b) moves, given results of a moves and b moves
    const auto prod{[N, inf{ranges::max(C | views::join) * K}](const auto& lhs, const auto& rhs) {
        // Initialize with a sufficiently large value
        vector ret(N, vector(N, inf));
        for (const auto i : views::iota(0U, N))
            for (const auto j : views::iota(0U, N))
                for (const auto k : views::iota(0U, N))
                    ret[i][k] = min(ret[i][k], lhs[i][j] + rhs[j][k]);
        return ret;
    }};

    auto ans(C);
    --K;

    // Fast exponentiation
    while (K) {
        if (K & 1)
            ans = prod(ans, C);
        C = prod(C, C);
        K /= 2;
    }

    for (const auto i : views::iota(0U, N))
        cout << ans[i][i] << endl;
    return 0;
}