POSTLLM - Editorial

PROBLEM LINK:

Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4

Author: sai1707 and weirdflexbutok
Editorialist: iceknight1093

DIFFICULTY:

Medium

PREREQUISITES:

Sorting, expected value, elementary combinatorics

PROBLEM:

There are N positions on the x-axis, represented by the array d.
You would like to choose M of them, after which you will choose a random permutation of these M points and move between them in this order.

Find the minimum expected distance you travel, along with the number of ways to choose points such that this minimum is achieved.

EXPLANATION:

Of course, we first try to compute the minimum.
For this, let’s understand how exactly the expected distance is computed.

Suppose we choose points \{x_1, x_2, \ldots, x_M\}, with x_i \leq x_{i+1}.
Since distance calculation is just the length of the straight path between points, the easiest way to calculate the expected distance travelled is using linearity of expectation: we can find, for each segment [x_i, x_{i+1}], the probability that this segment will be crossed.
If p_i is this probability, then the expected distance is simply \sum_{i=1}^{M-1} p_i \cdot (x_{i+1} - x_i).

So, let’s attempt to compute p_i.
[x_i, x_{i+1}] will be crossed if and only if the postman has to move from a point that’s \leq x_i to a point that’s \geq x_{i+1} (or vice versa).
There are 2\cdot i \cdot (M-i) ways to choose such a pair, after which this pair and the remaining (M-2) elements can be arranged in any of (M-1)! ways.
So, we obtain p_i = \frac{1}{M!} \cdot 2\cdot i \cdot (M-i) \cdot (M-1)! = \frac{2\cdot i\cdot (M-i)}{M} .

Our goal is to compute the answer multiplied by M, and M is a common denominator for all the p_i, so we can just discard the denominator entirely.
This leaves us with the value \sum_{i=1}^N 2i\cdot (M-i)\cdot (x_{i+1} - x_i), which we want to minimize.

Looking at this expression, it can be seen that to minimize this, it’s ideal for all the x_i to be “close together”.
In particular, the minimum value will always be obtained by some ‘window’ of d_i, i.e. by choosing the subset \{d_i, d_{i+1}, d_{i+2}, \ldots, d_{i+M-1}\} for some i - of course, after sorting the array d.


Our first task is hence to compute the minimum of \sum_{j=1}^M 2j\cdot (M-j) \cdot (d_{i+j} - d_{i+j-1}) across all valid i.

To do this, we’ll use a sort of sliding window.
Suppose this value is known for some i. We’ll try to recompute it for i+1 without explicitly doing \mathcal{O}(M) work.

To do this easily, observe what happens if we expand the summation out.
It can be observed that:

  • The coefficient of d_i is 2\cdot (1-M)
  • The coefficient of d_{i+1} is 2\cdot (3-M)
  • The coefficient of d_{i+2} is 2\cdot (5-M)
    \vdots
  • The coefficient of d_{i+M-1} is 2\cdot (M-1).

This allows us to update pretty easily when moving from i to i+1, because the only changes are:

  • The term d_i \cdot 2 \cdot (1 - M) disappears.
  • The term d_{i+M}\cdot 2 \cdot (M-1) appears.
  • For every other term, its coefficient reduces by exactly 2.
    So, the overall change in the sum is easily found by simply knowing the sum d_{i+1} + d_{i+2} + \ldots + d_{i+M-1}, which is just a range sum of the array d.

Since we can move from the window at i to the window at i+1 in constant time, we can now simply check all windows and compute the minimum value.


Once this minimum is known, let’s move to counting the number of configurations that achieve the minimum.

Let’s fix i, the starting point of the window.
If the window [i, i+M-1] doesn’t give us the minimum cost, we can ignore this i.
Otherwise, there’s at least one subset starting here with minimum cost.

Now, note that for optimality to hold, we really don’t have many options.
Suppose there are c_1 occurrences of d_i in this window, and c_2 occurrences of d_{i+M-1}.
Then, it doesn’t really matter which occurrences of them we take, so if their overall counts are n_1 and n_2 respectively, we have \binom{n_1}{c_1} \cdot \binom{n_2}{c_2} options to choose them.
In fact, this is all the freedom we have: all the ‘internal’ points must be included no matter what; otherwise the cost will increase.

So, the number of ways can be found by just adding up some products of binomial coefficients, which is easy enough to do.

Note that there is one edge case: when the minimum distance travelled is 0, we’ll have d_i = d_{i+M-1}, in which case rather than a product of binomial coefficients, the number of ways is just the number of ways of choosing M elements from all occurrences of d_i.
Make sure not to overcount here.

TIME COMPLEXITY:

\mathcal{O}(N\log N) per testcase.

CODE:

Author's code (C++)
#include<bits/stdc++.h>
using namespace std;
using i64 = long long;
// #define int long long
template <typename T>
T inverse(T a, T m) {
    T u = 0, v = 1;
    while (a != 0) {
        T t = m / a;
        m -= t * a; swap(a, m);
        u -= t * v; swap(u, v);
    }
    assert(m == 1);
    return u;
}
 
template <typename T>
class Modular {
public:
    using Type = typename decay<decltype(T::value)>::type;
 
    constexpr Modular() : value() {}
    template <typename U>
    Modular(const U& x) {
        value = normalize(x);
    }
 
    template <typename U>
    static Type normalize(const U& x) {
        Type v;
        if (-mod() <= x && x < mod()) v = static_cast<Type>(x);
        else v = static_cast<Type>(x % mod());
        if (v < 0) v += mod();
        return v;
    }
 
    const Type& operator()() const { return value; }
    template <typename U>
    explicit operator U() const { return static_cast<U>(value); }
    constexpr static Type mod() { return T::value; }
 
    Modular& operator+=(const Modular& other) { if ((value += other.value) >= mod()) value -= mod(); return *this; }
    Modular& operator-=(const Modular& other) { if ((value -= other.value) < 0) value += mod(); return *this; }
    template <typename U> Modular& operator+=(const U& other) { return *this += Modular(other); }
    template <typename U> Modular& operator-=(const U& other) { return *this -= Modular(other); }
    Modular& operator++() { return *this += 1; }
    Modular& operator--() { return *this -= 1; }
    Modular operator++(int) { Modular result(*this); *this += 1; return result; }
    Modular operator--(int) { Modular result(*this); *this -= 1; return result; }
    Modular operator-() const { return Modular(-value); }
 
    template <typename U = T>
    typename enable_if<is_same<typename Modular<U>::Type, int>::value, Modular>::type & operator*=(const Modular& rhs) {
        value = normalize(static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value));
        return *this;
    }
    template <typename U = T>
    typename enable_if<is_same<typename Modular<U>::Type, long long>::value, Modular>::type & operator*=(const Modular& rhs) {
        long long q = static_cast<long long>(static_cast<long double>(value) * rhs.value / mod());
        value = normalize(value * rhs.value - q * mod());
        return *this;
    }
    template <typename U = T>
    typename enable_if < !is_integral<typename Modular<U>::Type>::value, Modular >::type & operator*=(const Modular& rhs) {
        value = normalize(value * rhs.value);
        return *this;
    }
 
    Modular& operator/=(const Modular& other) { return *this *= Modular(inverse(other.value, mod())); }
 
    friend const Type& abs(const Modular& x) { return x.value; }
 
    template <typename U>
    friend bool operator==(const Modular<U>& lhs, const Modular<U>& rhs);
 
    template <typename U>
    friend bool operator<(const Modular<U>& lhs, const Modular<U>& rhs);
 
    template <typename V, typename U>
    friend V& operator>>(V& stream, Modular<U>& number);
 
private:
    Type value;
};
 
template <typename T> bool operator==(const Modular<T>& lhs, const Modular<T>& rhs) { return lhs.value == rhs.value; }
template <typename T, typename U> bool operator==(const Modular<T>& lhs, U rhs) { return lhs == Modular<T>(rhs); }
template <typename T, typename U> bool operator==(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) == rhs; }
 
template <typename T> bool operator!=(const Modular<T>& lhs, const Modular<T>& rhs) { return !(lhs == rhs); }
template <typename T, typename U> bool operator!=(const Modular<T>& lhs, U rhs) { return !(lhs == rhs); }
template <typename T, typename U> bool operator!=(U lhs, const Modular<T>& rhs) { return !(lhs == rhs); }
 
template <typename T> bool operator<(const Modular<T>& lhs, const Modular<T>& rhs) { return lhs.value < rhs.value; }
 
template <typename T> Modular<T> operator+(const Modular<T>& lhs, const Modular<T>& rhs) { return Modular<T>(lhs) += rhs; }
template <typename T, typename U> Modular<T> operator+(const Modular<T>& lhs, U rhs) { return Modular<T>(lhs) += rhs; }
template <typename T, typename U> Modular<T> operator+(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) += rhs; }
 
template <typename T> Modular<T> operator-(const Modular<T>& lhs, const Modular<T>& rhs) { return Modular<T>(lhs) -= rhs; }
template <typename T, typename U> Modular<T> operator-(const Modular<T>& lhs, U rhs) { return Modular<T>(lhs) -= rhs; }
template <typename T, typename U> Modular<T> operator-(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) -= rhs; }
 
template <typename T> Modular<T> operator*(const Modular<T>& lhs, const Modular<T>& rhs) { return Modular<T>(lhs) *= rhs; }
template <typename T, typename U> Modular<T> operator*(const Modular<T>& lhs, U rhs) { return Modular<T>(lhs) *= rhs; }
template <typename T, typename U> Modular<T> operator*(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) *= rhs; }
 
template <typename T> Modular<T> operator/(const Modular<T>& lhs, const Modular<T>& rhs) { return Modular<T>(lhs) /= rhs; }
template <typename T, typename U> Modular<T> operator/(const Modular<T>& lhs, U rhs) { return Modular<T>(lhs) /= rhs; }
template <typename T, typename U> Modular<T> operator/(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) /= rhs; }
 
template<typename T, typename U>
Modular<T> power(const Modular<T>& a, const U& b) {
    assert(b >= 0);
    Modular<T> x = a, res = 1;
    U p = b;
    while (p > 0) {
        if (p & 1) res *= x;
        x *= x;
        p >>= 1;
    }
    return res;
}
 
template <typename T>
bool IsZero(const Modular<T>& number) {
    return number() == 0;
}
 
template <typename T>
string to_string(const Modular<T>& number) {
    return to_string(number());
}
 
// U == std::ostream? but done this way because of fastoutput
template <typename U, typename T>
U& operator<<(U& stream, const Modular<T>& number) {
    return stream << number();
}
 
// U == std::istream? but done this way because of fastinput
template <typename U, typename T>
U& operator>>(U& stream, Modular<T>& number) {
    typename common_type<typename Modular<T>::Type, long long>::type x;
    stream >> x;
    number.value = Modular<T>::normalize(x);
    return stream;
}
 
/*
using ModType = int;
 
struct VarMod { static ModType value; };
ModType VarMod::value;
ModType& md = VarMod::value;
using Mint = Modular<VarMod>;
*/
 
constexpr int md = 998244353;
using Mint = Modular<std::integral_constant<decay<decltype(md)>::type, md>>;
 
vector<Mint> fact(1, 1);
vector<Mint> inv_fact(1, 1);
 
Mint C(int n, int k) {
    if (k < 0 || k > n) {
        return 0;
    }
    while ((int) fact.size() < n + 1) {
        fact.push_back(fact.back() * (int) fact.size());
        inv_fact.push_back(1 / fact.back());
    }
    return fact[n] * inv_fact[k] * inv_fact[n - k];
}
int32_t main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(0);
    int T;
    cin >> T;
    while (T--) {
        int n, m;
        cin >> m >> n;
        vector<i64> d(m);
        vector<i64> pre(m + 1);
        for (int i = 0; i < m; i++) {
            cin >> d[i];
        }
        sort(d.begin(), d.end());
        i64 sm = 0, sm1 = 0;
        int factor = n + 1;
        map<i64, i64> cnt;
        for (int i = 0; i < m; i++) {
            pre[i + 1] = pre[i] + d[i];
            cnt[d[i]]++;
        }
        for (int i = n - 1; i >= 0; i--) {
            factor -= 2;
            sm += factor * d[i];
        }
        factor = -factor;
        i64 ans = sm;
 
        int l = std::upper_bound(d.begin(), d.end(), d[0]) - d.begin() - 1;
        int r = std::upper_bound(d.begin(), d.end(), d[n - 1] - 1) - d.begin();
        int lc = l - (n - 1 - n + 1) + 1;
        int rc = (n - 1) - r + 1;
 
        array<i64, 4> info = {lc, rc, d[0], d[n - 1]};
        map<array<i64, 4>, int> anss;
        anss[info]++;
 
        for (int i = n; i < m; i++) {
            sm -= 2 * (pre[i] - pre[i - n]);
            sm += (n - 1) * d[i];
            sm += (factor + 2) * d[i - n];
 
            int l = std::upper_bound(d.begin(), d.end(), d[i - n + 1]) - d.begin() - 1;
            int r = std::upper_bound(d.begin(), d.end(), d[i] - 1) - d.begin();
 
            int lc = l - (i - n + 1) + 1;
            int rc = i - r + 1;
 
            array<i64, 4> info = {lc, rc, d[i - n + 1], d[i]};
 
            if (sm == ans) {
                anss[info]++;
            }
            if (sm < ans) {
                anss.clear();
                anss[info]++;
            }
            ans = std::min(ans, sm);
        }
        Mint ways = 0;
        set<int> done;
        cout << ans * 2 << ' ';
        for (auto &[x, y] : anss) {
            auto [lc, rc, a, b] = x;
            if (a == b and !done.count(a)) {
                ways += C(cnt[a], n);
                done.insert(a);
            }
            else if (a != b) {
                ways += C(cnt[a], lc) * C(cnt[b], rc);
            }
        }
        cout << ways << '\n';
    }
}

Editorialist's code (C++)
// #include <bits/allocator.h>
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include "bits/stdc++.h"
using namespace std;
using ll = long long int;
mt19937_64 RNG(chrono::high_resolution_clock::now().time_since_epoch().count());
 
/**
 * Integers modulo p, where p is a prime
 * Source: Aeren (modified from tourist?)
 *         Modmul for 64-bit mod from kactl:ModMulLL
 * Works with p < 7.2e18 with x87 80-bit long double, and p < 2^52 ~ 4.5e12 with 64-bit
 */
template<typename T>
struct Z_p{
	using Type = typename decay<decltype(T::value)>::type;
	static vector<Type> MOD_INV;
	constexpr Z_p(): value(){ }
	template<typename U> Z_p(const U &x){ value = normalize(x); }
	template<typename U> static Type normalize(const U &x){
		Type v;
		if(-mod() <= x && x < mod()) v = static_cast<Type>(x);
		else v = static_cast<Type>(x % mod());
		if(v < 0) v += mod();
		return v;
	}
	const Type& operator()() const{ return value; }
	template<typename U> explicit operator U() const{ return static_cast<U>(value); }
	constexpr static Type mod(){ return T::value; }
	Z_p &operator+=(const Z_p &otr){ if((value += otr.value) >= mod()) value -= mod(); return *this; }
	Z_p &operator-=(const Z_p &otr){ if((value -= otr.value) < 0) value += mod(); return *this; }
	template<typename U> Z_p &operator+=(const U &otr){ return *this += Z_p(otr); }
	template<typename U> Z_p &operator-=(const U &otr){ return *this -= Z_p(otr); }
	Z_p &operator++(){ return *this += 1; }
	Z_p &operator--(){ return *this -= 1; }
	Z_p operator++(int){ Z_p result(*this); *this += 1; return result; }
	Z_p operator--(int){ Z_p result(*this); *this -= 1; return result; }
	Z_p operator-() const{ return Z_p(-value); }
	template<typename U = T>
	typename enable_if<is_same<typename Z_p<U>::Type, int>::value, Z_p>::type &operator*=(const Z_p& rhs){
		#ifdef _WIN32
		uint64_t x = static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value);
		uint32_t xh = static_cast<uint32_t>(x >> 32), xl = static_cast<uint32_t>(x), d, m;
		asm(
			"divl %4; \n\t"
			: "=a" (d), "=d" (m)
			: "d" (xh), "a" (xl), "r" (mod())
		);
		value = m;
		#else
		value = normalize(static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value));
		#endif
		return *this;
	}
	template<typename U = T>
	typename enable_if<is_same<typename Z_p<U>::Type, int64_t>::value, Z_p>::type &operator*=(const Z_p &rhs){
		uint64_t ret = static_cast<uint64_t>(value) * static_cast<uint64_t>(rhs.value) - static_cast<uint64_t>(mod()) * static_cast<uint64_t>(1.L / static_cast<uint64_t>(mod()) * static_cast<uint64_t>(value) * static_cast<uint64_t>(rhs.value));
		value = normalize(static_cast<int64_t>(ret + static_cast<uint64_t>(mod()) * (ret < 0) - static_cast<uint64_t>(mod()) * (ret >= static_cast<uint64_t>(mod()))));
		return *this;
	}
	template<typename U = T>
	typename enable_if<!is_integral<typename Z_p<U>::Type>::value, Z_p>::type &operator*=(const Z_p &rhs){
		value = normalize(value * rhs.value);
		return *this;
	}
	template<typename U>
	Z_p &operator^=(U e){
		if(e < 0) *this = 1 / *this, e = -e;
		Z_p res = 1;
		for(; e; *this *= *this, e >>= 1) if(e & 1) res *= *this;
		return *this = res;
	}
	template<typename U>
	Z_p operator^(U e) const{
		return Z_p(*this) ^= e;
	}
	Z_p &operator/=(const Z_p &otr){
		Type a = otr.value, m = mod(), u = 0, v = 1;
		if(a < (int)MOD_INV.size()) return *this *= MOD_INV[a];
		while(a){
			Type t = m / a;
			m -= t * a; swap(a, m);
			u -= t * v; swap(u, v);
		}
		assert(m == 1);
		return *this *= u;
	}
	template<typename U> friend const Z_p<U> &abs(const Z_p<U> &v){ return v; }
	Type value;
};
template<typename T> bool operator==(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value == rhs.value; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator==(const Z_p<T>& lhs, U rhs){ return lhs == Z_p<T>(rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator==(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) == rhs; }
template<typename T> bool operator!=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return !(lhs == rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator!=(const Z_p<T> &lhs, U rhs){ return !(lhs == rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator!=(U lhs, const Z_p<T> &rhs){ return !(lhs == rhs); }
template<typename T> bool operator<(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value < rhs.value; }
template<typename T> bool operator>(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value > rhs.value; }
template<typename T> bool operator<=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value <= rhs.value; }
template<typename T> bool operator>=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value >= rhs.value; }
template<typename T> Z_p<T> operator+(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator+(const Z_p<T> &lhs, U rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator+(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T> Z_p<T> operator-(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator-(const Z_p<T>& lhs, U rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator-(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T> Z_p<T> operator*(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator*(const Z_p<T>& lhs, U rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator*(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T> Z_p<T> operator/(const Z_p<T> &lhs, const Z_p<T> &rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator/(const Z_p<T>& lhs, U rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator/(U lhs, const Z_p<T> &rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T> istream &operator>>(istream &in, Z_p<T> &number){
	typename common_type<typename Z_p<T>::Type, int64_t>::type x;
	in >> x;
	number.value = Z_p<T>::normalize(x);
	return in;
}
template<typename T> ostream &operator<<(ostream &out, const Z_p<T> &number){ return out << number(); }
 
/*
using ModType = int;
struct VarMod{ static ModType value; };
ModType VarMod::value;
ModType &mod = VarMod::value;
using Zp = Z_p<VarMod>;
*/
 
// constexpr int mod = 1e9 + 7; // 1000000007
constexpr int mod = (119 << 23) + 1; // 998244353
// constexpr int mod = 1e9 + 9; // 1000000009
using Zp = Z_p<integral_constant<decay<decltype(mod)>::type, mod>>;
 
template<typename T> vector<typename Z_p<T>::Type> Z_p<T>::MOD_INV;
template<typename T = integral_constant<decay<decltype(mod)>::type, mod>>
void precalc_inverse(int SZ){
	auto &inv = Z_p<T>::MOD_INV;
	if(inv.empty()) inv.assign(2, 1);
	for(; inv.size() <= SZ; ) inv.push_back((mod - 1LL * mod / (int)inv.size() * inv[mod % (int)inv.size()]) % mod);
}
 
template<typename T>
vector<T> precalc_power(T base, int SZ){
	vector<T> res(SZ + 1, 1);
	for(auto i = 1; i <= SZ; ++ i) res[i] = res[i - 1] * base;
	return res;
}
 
template<typename T>
vector<T> precalc_factorial(int SZ){
	vector<T> res(SZ + 1, 1); res[0] = 1;
	for(auto i = 1; i <= SZ; ++ i) res[i] = res[i - 1] * i;
	return res;
}
 
int main()
{
    ios::sync_with_stdio(false); cin.tie(0);

    auto fac = precalc_factorial<Zp>(1000006);
    auto inv = fac;
    inv.back() = Zp(1) / inv.back();
    for (int i = inv.size() - 2; i >= 0; --i) inv[i] = (i+1) * inv[i+1];
    auto C = [&] (int n, int r) {
        if (n < r or r < 0) return Zp(0);
        return fac[n] * inv[r] * inv[n-r];
    };
 
    int t; cin >> t;
    while (t--) {
        int n, m; cin >> n >> m;
        vector<int> d(n);
        map<int, int> freq;
        for (int &x : d) {
            cin >> x;
            ++freq[x];
        }

        ranges::sort(d);
        ll ans = 0;
        ll cur = 0, sm = 0;
        for (int i = 0; i < m; ++i) {
            sm += d[i];
            cur += 1ll*d[i]*(-m+1+2*i);
        }
        ans = cur;
        Zp ways = 0;

        {
            auto it = ranges::lower_bound(d, d[m-1]) - begin(d);
            int ct = m - it;
            ways = C(freq[d[m-1]], ct);
        }
        
        for (int i = m; i < n; ++i) {
            sm -= d[i-m];
            cur -= 1ll*d[i-m]*(-m+1);
            cur -= 2*sm;
            cur += 1ll*d[i]*(m-1);
            sm += d[i];

            if (ans < cur) continue;
            if (ans > cur) {
                ans = cur;
                ways = 0;
            }

            int left = (ranges::upper_bound(d, d[i-m+1]) - begin(d)) - (i-m+1);
            int right = i+1 - (ranges::lower_bound(d, d[i]) - begin(d));
            ways += C(freq[d[i-m+1]], left) * C(freq[d[i]], right);
        }
        if (ans == 0) {
            ways = 0;
            for (auto [x, y] : freq) {
                ways += C(y, m);
            }
        }
        cout << 2*ans << ' ' << ways << '\n';
    }
}

Solved this one within 5 minutes after the end. What a pity!