详细请查看注释

参考代码：

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const int N = 25, INF = 0x3f3f3f3f;

int n, m, R[N], H[N];
// 分别存储顶层到当前层体积最小值和面积最小值（包括当前层）
int minv[N], mins[N], res = INF;

/*
    注意：无需考虑π的存在！

    剪枝：
        1. 优化搜索书顺序：自底向上枚举，自大到小枚举R, H, 先枚举R再枚举H（R计算时是平方级别）
        2. 可行性剪枝：（记顶层为第一层）

            设底层到当前第u层（不包括第u层）总体积为v，总体积为n：则 n - v >= R[u]^2 * H[u]
                当H[u]取1时R[u]达到最大，R[u] <= sqrt(n - v)
                                         H[u] <= (n - v / R[u]^2)

            R[u]范围：u <= R[u] <= min(R[u + 1] - 1, sqrt(n - v))
            H[u]范围：u <= H[u] <= min(H[u + 1] - 1, (n - v) / R[u]^2)

        3. 可行性剪枝：
            minv[u], mins[u] 分别存储顶层到当前层体积最小值和面积最小值（包括当前层）
            res表示当前表面积总和的最优解！s表示当前层的表面积！

            1. minv[u] + v <= n
            2. mins[u] + s < res  ---> 最优性剪枝

        4. 通过表面积体积公式推导一个关系：
            S[1-u] = 2R[k] * H[k] (k = 1,2,3...u)
            n - v = R[k]^2 * H[k] (k = 1,2,3...u)

            给s[1-u]除一次R[u + 1]乘一次R[u + 1]：（R[u + 1] >= R[k]）
                s[1-u] = (2 / R[u + 1]) * R[k] * H[k] * R[u + 1] 
                       >= (2 / R[u + 1]) * R[k]^2 * H[k]
                       >= (2 / R[u + 1]) * (n - v)

            因此：s + s[1-u] >= res 
                即s + 2(n - v) / R[u + 1] >= res  ---> 直接return
*/

// 当前层编号 当前总体积 当前总表面积
void dfs(int u, int v, int s){
    if(minv[u] + v > n) return;
    if(mins[u] + s >= res) return;
    if(s + 2 * (n - v) / R[u + 1] >= res) return;

    if(!u) {
        if(v == n) res = s;
        return;
    }

    // 先枚举r在枚举h, 合法范围内枚举
    for(int r = min(R[u + 1] - 1, (int)sqrt(n - v)); r >= u; r --)
        for(int h = min(H[u + 1] - 1, (n - v) / r / r); h >= u; h --){
            int t = 0;
            // 若是底层额外加上底层底面积（即每个顶层顶面积的和）
            if(u == m) t = r * r;
            R[u] = r, H[u] = h;
            dfs(u - 1, v + r * r * h, s + 2 * r * h + t);
        }
}

int main(){
    cin >> n >> m;
    // 预处理 顶层到当前层体积最小值和面积最小值（包括当前层）
    for(int i = 1; i <= m; i ++){
        minv[i] = minv[i - 1] + i * i * i;
        mins[i] = mins[i - 1] + 2 * i * i;
    }

    R[m + 1] = H[m + 1] = INF;
    dfs(m, 0, 0);
    if(res == INF) cout << 0 << endl;
    else cout << res << endl;
    return 0;
}