※ 参加記や問題そのものの解法ではありません

頑丈さがわかっているとき、現実的な時間で厳密解を求めることができます。

掘削のパワーを \(S_{i,j}\) とすればよいのは明らかです。よって、どこを掘削するかが問題になります。

まず、グリッドを \(N^2\) 頂点のグラフと捉えます。ここで、隣接する頂点間に辺があり、頂点 \((i, j)\) の重みは \(c_{(i,j)} = C+S_{i,j}\) であるとします。さらに重み \(0\) の頂点 \(w\) を加えて、\(w\) とすべての水源 \((a_i,b_i)\) の間に辺を張ります。

求めるものは、\(w\) とすべての家 \((c_i, d_i)\) を含む部分グラフであって頂点の重みの和が最小になるものです。

ここで、\(S\) を \(\{(c_0,d_0),\ldots,(c_{K-1},d_{K-1})\}\) の部分集合として、次の動的計画法を考えます。

\(\mathrm{dp}[S][v]\)：頂点集合に \(S\cup\{v\}\) を含む部分グラフの \(v\) 以外の頂点の重みの和の最小値

明らかに \(\mathrm{dp}[\{(c_i, d_i)\}][(c_i, d_i)] = 0\) であり、それ以外は以下の漸化式に従って求めることができます。

\(\mathrm{dp}[S][v]=\min(\{\mathrm{dp}[S][u]+c_u\mid (u,v)\text{は辺} \} \\\phantom{\mathrm{dp}[S][v]=\min(}{}\cup\{\mathrm{dp}[T][v]+\mathrm{dp}[S\setminus T][v]\mid T\subsetneq S\})\)

実際には、\(S\) ごとに \(\mathrm{dp}[T][v]+\mathrm{dp}[S\setminus T][v]\) による遷移を終えてから全頂点を始点とした Dijkstra 法で \(\mathrm{dp}[S][u]+c_u\) による遷移をすればよいです。

また、\(\mathrm{dp}[S][v]\) を実現する方法を得るため、\(v\) と接続する頂点を \(S\) や \(T\) や \(S\setminus T\) とともに記録しておく必要があります。

時間計算量は \(O(3^KN^2+2^KN^2\log N^2)=O((3^K+2^K\log N)N^2)\) です。

seed = 0 について、手元では約 30 秒で以下の解（Cost = 67275）を得ることができました。

実装例（Rust；テストケースをそのまま入力する）

use itertools::{chain, iproduct, izip};
use proconio::input;
use std::{cmp::Reverse, collections::BinaryHeap};

fn main() {
    input!(
        n: usize,
        n_water: usize,
        n_house: usize,
        cost0: u32,
        mut cost: [[u32; n]; n],
        water: [(usize, usize); n_water],
        house: [(usize, usize); n_house]
    );

    cost.iter_mut().flatten().for_each(|c| *c += cost0);

    let mut dp = vec![vec![vec![std::u32::MAX; n]; n]; 1 << n_house];
    let mut prev = vec![vec![vec![vec![]; n]; n]; 1 << n_house];

    let mut dp_w = vec![std::u32::MAX; 1 << n_house];
    let mut prev_w = vec![vec![]; 1 << n_house];

    for (k, &(i, j)) in izip!(0.., &house) {
        dp[1 << k][i][j] = 0;
    }

    for set in 1usize..1 << n_house {
        let mut subset = set.wrapping_sub(1) & set;
        while subset != 0 {
            for (i, j) in iproduct!(0..n, 0..n) {
                let dp_ = dp[subset][i][j] + dp[set & !subset][i][j];
                if dp[set][i][j] > dp_ {
                    dp[set][i][j] = dp_;
                    prev[set][i][j] = chain(&prev[subset][i][j], &prev[set & !subset][i][j])
                        .copied()
                        .collect();
                }
            }
            let dp_w_ = dp_w[subset] + dp_w[set & !subset];
            if dp_w[set] > dp_w_ {
                dp_w[set] = dp_w_;
                prev_w[set] = chain(&prev_w[subset], &prev_w[set & !subset])
                    .copied()
                    .collect();
            }
            subset = subset.wrapping_sub(1) & set;
        }
        let dp = &mut dp[set];
        let prev = &mut prev[set];
        let mut heap = iproduct!(0..n, 0..n)
            .map(|(i, j)| (Reverse(dp[i][j]), (i, j)))
            .collect::<BinaryHeap<_>>();
        while let Some((Reverse(d1), (i1, j1))) = heap.pop() {
            if dp[i1][j1] != d1 {
                continue;
            }
            for (i2, j2) in adj4((i1, j1), n) {
                let d2_ = dp[i1][j1] + cost[i1][j1];
                if dp[i2][j2] > d2_ {
                    dp[i2][j2] = d2_;
                    prev[i2][j2] = vec![(set, (i1, j1))];
                    heap.push((Reverse(d2_), (i2, j2)));
                }
            }
        }
        for &(i, j) in &water {
            let dp_w_ = dp[i][j] + cost[i][j];
            if dp_w[set] > dp_w_ {
                dp_w[set] = dp_w_;
                prev_w[set] = vec![(set, (i, j))];
            }
        }
    }

    let all = !(!0 << n_house);

    let mut stack = prev_w[all].clone();
    while let Some((set, (i, j))) = stack.pop() {
        println!("{} {} {}", i, j, cost[i][j]);
        stack.extend(&prev[set][i][j]);
    }
}

const N1: usize = 1usize.wrapping_neg();
const D4: [(usize, usize); 4] = [(1, 0), (0, 1), (N1, 0), (0, N1)];

fn adj((i, j): (usize, usize), dir: usize) -> (usize, usize) {
    let (di, dj) = D4[dir];
    (i.wrapping_add(di), j.wrapping_add(dj))
}

fn adj4(ix: (usize, usize), n: usize) -> impl Iterator<Item = (usize, usize)> {
    (0..4)
        .map(move |dir| adj(ix, dir))
        .filter(move |&(i, j)| i < n && j < n)
}