妖精の演奏 [yukicoder No.572] - はまやんはまやんはまやん

前提知識

ダブリング

解法

非想定解（落とされるかも）
https://yukicoder.me/submissions/208104
キレイじゃない解法。
　
始点sと終点tと中継点cを決めて、その間の距離を全探索する。
最適な解は「s -> ... -> c -> ... -> c -> ... -> c -> ... -> t」の様な遷移になると仮定した解法。
これは最初と最後以外は効率のいい所をぐるぐるまわるだろうという考察からなる。
最初は「rep(s, 0, M) rep(c, 0, M) rep(t, 0, M) rep(len1, 1, 32) rep(len2, 3, 32) rep(len3, 1, 32)」というループでやっていた。
これは、点sから長さlen1で点cまで行き、点cから点cまで長さlen2のループが任意回あり、点cから点tまで長さlen3で進んで終わりというパターンを全探索している。
途中のループがちゃんとできるかは、(N - len1 - len3 + 1)が(len2 - 1)で割り切れるかで判定している。
　
あとは、途中の長さを予めDPで計算しておけば、計算量は最大2^30(=10^9)くらい。
dp[len][s][t] := 点sから初めて点tまで長さlenでの最大コスト
ギリギリ間に合いそうだが、最初は間に合わなかったので、len3は全てループではなく、ありえるものだけにすると大分早くなって間に合った。

typedef long long ll;
ll N; int M;
// dp[len][s][t] := max
ll dp[1010][30][30];
ll A[30][30];
//---------------------------------------------------------------------------------------------------
void _main() {
    cin >> N >> M;
    rep(y, 0, M) rep(x, 0, M) cin >> A[y][x];

    rep(s, 0, M) rep(t, 0, M) dp[2][s][t] = A[s][t];
    rep(le, 2, 1000) rep(s, 0, M) rep(t, 0, M) {
        rep(to, 0, M) dp[le + 1][s][to] = max(dp[le + 1][s][to], dp[le][s][t] + A[t][to]);
    }

    ll ans = 0;
    if (N <= 100) {
        rep(s, 0, M) rep(t, 0, M) ans = max(ans, dp[N][s][t]);
    } else {
        rep(s, 0, M) rep(c, 0, M) rep(t, 0, M) rep(len1, 1, 32) rep(len2, 3, 32) {
            if (len1 == 1 && s != c) continue;
            ll _d = N - len1 + 1;
            int m = _d % (len2 - 1);
            for (int len3 = len2 - 1 - m; len3 < 32; len3 += len2 - 1) {
                if (len3 == 1 && c != t) continue;
                ll d = N - len1 - len3 + 1;
                ll cc = dp[len1][s][c] + dp[len3][c][t];
                //assert(d % (len2 - 1) == 0);
                cc += dp[len2][c][c] * (d / (len2 - 1));
                if (ans < cc) {
                    ans = cc;
                }
            }
        }
    }

    cout << ans << endl;
}

　
想定解
https://yukicoder.me/submissions/208148

想定解はダブリング。
dp[d][x][y] := 長さ2^dで点xから点yまでの最大コスト
これはdp[d][x][y] = max{z=0..M-1}(dp[d-1][x][z]+dp[d-1][z][y])で計算できる。
あとは、(N-1)回遷移するときに、bitが立っていればdpを使って遷移をしていけばよい。
(N-1なのは一番最初の点を決めるがこれは遷移回数では無いため引く)
　
ans2[x][y] := 現時点で点xから点yへ遷移するときの最大コスト
これを更新していくが、例えばN-1=5であれば、5=1+4であるため、2^0の遷移と2^2の遷移を組み合わせればよい。
以下のコードを見てもらったほうが何をしているか分かりやすいかもしれない。

typedef long long ll;
ll N; int M;
ll A[30][30];
ll dp[40][30][30];
ll ans[30][30];
//---------------------------------------------------------------------------------------------------
void _main() {
    cin >> N >> M;
    rep(y, 0, M) rep(x, 0, M) cin >> A[y][x];

    rep(y, 0, M) rep(x, 0, M) dp[0][y][x] = A[y][x];
    rep(d, 1, 40) rep(x, 0, M) rep(y, 0, M) {
        rep(z, 0, M) dp[d][x][y] = max(dp[d][x][y], dp[d - 1][x][z] + dp[d - 1][z][y]);
    }
    N--;
    int d = 0;
    while (N) {
        if (N % 2) {
            ll ans2[30][30];
            rep(x, 0, M) rep(y, 0, M) ans2[x][y] = 0;

            rep(x, 0, M) rep(y, 0, M) {
                rep(z, 0, M) ans2[x][y] = max(ans2[x][y], ans[x][z] + dp[d][z][y]);
            }

            swap(ans, ans2);
        }

        N /= 2;
        d++;
    }

    ll a = 0;
    rep(x, 0, M) rep(y, 0, M) a = max(a, ans[x][y]);
    cout << a << endl;
}