\(A_1,A_2,\dots,A_N\) から重複を除いて，現れる値が \(y_1 < y_2 < \dots < y_L\) になったとします． \(S_j\) を \(A_i\in(0,y_j]\) となる \(i\) の個数とします．

\[dp[i][j] := (0,y_j] を i 個の区間に分けたときの最小値\]

とした動的計画法を考えます．遷移は

\[dp[i][j]=\min_{k<j} \left\lbrace dp[i-1][k]+y_j\times(S_j-S_k) \right\rbrace\]

となります．これは \(i\) ごとに convex hull trick を用いて \(O(NK)\) 時間となります．

ところで，コスト関数 \(w(k,j)=y_j\times(S_j-S_k)\) は Monge 性を満たします．

よって，これは Monge グラフ上の \(d\) -辺最短路長を計算するアルゴリズム (いわゆる Alien DP) で \(O(N \log (N \max A))\) 時間で解けます．

実装例

上では説明の都合上 \(A\) の重複を除きましたが，実装時には重複を除く必要はありません．

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

template <typename T>
class ConvexHullTrick {
   public:
    void add(T a, T b) {
        Line line(a, b);
        while (lines.size() >= 2 &&
               check(*(lines.end() - 2), lines.back(), line)) {
            lines.pop_back();
        }
        lines.push_back(line);
    }

    T get(T x) {
        while (lines.size() >= 2 && lines.front()(x) > lines[1](x)) {
            lines.pop_front();
        }
        return lines.front()(x);
    }

   private:
    struct Line {
        T a, b;
        Line(T a, T b) : a(a), b(b) {}
        T operator()(T x) const { return a * x + b; }
    };

    std::deque<Line> lines;

    static bool check(Line l1, Line l2, Line l3) {
        if (l2.a == l3.a) return l2.b >= l3.b;
        return 1.0 * (l2.b - l1.b) / (l2.a - l1.a) <=
               1.0 * (l3.b - l2.b) / (l3.a - l2.a);
    }
};

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    cout << fixed << setprecision(15);

    int N, K;
    cin >> N >> K;
    K = min(K, N);
    vector<ll> A(N);
    for (auto& x : A) cin >> x;
    sort(A.begin(), A.end());

    auto calc = [&](ll lambda) {
        vector<ll> dp(N + 1);
        ConvexHullTrick<ll> cht;
        cht.add(0, 0);
        for (int j = 1; j <= N; ++j) {
            dp[j] = cht.get(A[j - 1]) + A[j - 1] * j + lambda;
            cht.add(-j, dp[j]);
        }
        return dp[N] - lambda * K;
    };

    ll lb = -3e10, ub = 3e10;
    while (ub - lb > 2) {
        ll m1 = (2 * lb + ub) / 3;
        ll m2 = (lb + 2 * ub) / 3;
        ll v1 = calc(m1);
        ll v2 = calc(m2);
        if (v1 > v2) {
            ub = m2;
        } else {
            lb = m1;
        }
    }
    cout << max(calc(lb), calc(lb + 1)) << endl;
}