Monge 行列の min-plus 合成 (in O(N^2))
(Optimization/monge_min_plus_composition.hpp)

• View this file on GitHub
• Last update: 2026-11-12 11:56:29+09:00
• Include: #include "Optimization/monge_min_plus_composition.hpp"

概要

上三角 $(N+1) \times (N+1)$ Monge 行列 $A, B$ に対して

C[i][j] = min_{i <= k < j} (A[i][k] + A[k][j])

を高速に実行する。

いわゆる「斜めに埋めるテク」を用いることで、$O(N^2)$ の計算量で実行できる。

問題例

参考文献

• tatyam: Monge の手引き書
• 週刊 spaghetti_source: Knuth-Yao speedup

Required by

• Monge 全部乗せ (Optimization/monge_all.hpp)

Code

#pragma once
#include <bits/stdc++.h>
#include <cassert>
using namespace std;

// C[i][j] = min_{i <= k < j} (A[i][k] + B[k][j])
// A, B: 上三角 (N+1) x (N+1) Monge
template<class VAL, class FUNCA, class FUNCB> vector<vector<VAL>> min_plus_composition
(int N, const FUNCA &A, const FUNCB &B) {
    VAL INF = numeric_limits<VAL>::max() / 2;
    vector<vector<VAL>> res(N+1, vector<VAL>(N+1, INF));
    vector<int> argmin(N+1);
    for (int i = 0; i <= N; i++) res[i][i] = A(i, i) + B(i, i), argmin[i] = i;
    for (int s = 1; s <= N; s++) {
        vector<int> argmin2(N + 1 - s);
        for (int i = 0; i <= N - s; i++) {
            int j = i + s, p = argmin[i], q = argmin[i + 1];
            for (int k = p; k <= q; k++) {
                if (res[i][j] > A(i, k) + B(k, j)) {
                    res[i][j] = A(i, k) + B(k, j);
                    argmin2[i] = k;
                }
            }
        }
        swap(argmin, argmin2);
    }
    return res;
}

#line 2 "Optimization/monge_min_plus_composition.hpp"
#include <bits/stdc++.h>
#line 4 "Optimization/monge_min_plus_composition.hpp"
using namespace std;

// C[i][j] = min_{i <= k < j} (A[i][k] + B[k][j])
// A, B: 上三角 (N+1) x (N+1) Monge
template<class VAL, class FUNCA, class FUNCB> vector<vector<VAL>> min_plus_composition
(int N, const FUNCA &A, const FUNCB &B) {
    VAL INF = numeric_limits<VAL>::max() / 2;
    vector<vector<VAL>> res(N+1, vector<VAL>(N+1, INF));
    vector<int> argmin(N+1);
    for (int i = 0; i <= N; i++) res[i][i] = A(i, i) + B(i, i), argmin[i] = i;
    for (int s = 1; s <= N; s++) {
        vector<int> argmin2(N + 1 - s);
        for (int i = 0; i <= N - s; i++) {
            int j = i + s, p = argmin[i], q = argmin[i + 1];
            for (int k = p; k <= q; k++) {
                if (res[i][j] > A(i, k) + B(k, j)) {
                    res[i][j] = A(i, k) + B(k, j);
                    argmin2[i] = k;
                }
            }
        }
        swap(argmin, argmin2);
    }
    return res;
}

Back to top page