Let \(M\) denote the maximum value for the input coordinates. We perform a binary search for the sought value, solving the problem “Can we draw a square with side length \(S\) that contains all the points?” \(O(\log M)\) times.

Focusing only on the \(x\) coordinates, we get the following problem:

(Problem) There are \(N\) points at coordinates \(x_1, \dots, x_N\) on a number line. Choose a closed interval \([l, l+S]\) and move all \(N\) points inside it (the same operations as in the original problem are available). Find the minimum total cost required for this.

Solve this problem separately for the \(x\)-coordinates and the \(y\)-coordinates, and if the total answer is at most \(K\), then the answer to the initial problem is Yes; otherwise, it is No.

All that remains is to choose \(l\) appropriately. The optimal choice is the \(N\)-th value when the \(2N\) values \(x_1, x_2, \dots, x_N, x_1 - S, x_2 - S, \dots, x_N - S\) are sorted in ascending order (precisely, it can be any value between the \(N\)-th and the \((N+1)\)-th values). This can be demonstrated by considering the changes in the slope of the piecewise linear graph obtained by plotting the total cost against \(l\) (the slope is \(-N\) for sufficiently small \(l\), and as \(l\) increases, it increases by \(1\) at each of the \(2N\) values, and the points corresponding to slope \(0\) are optimal).

If you sort \(x_1, x_2, \dots, x_N\) beforehand, “the \(N\)-th value when the \(2N\) values \(x_1, x_2, \dots, x_N, x_1 - S, x_2 - S, \dots, x_N - S\) are sorted in ascending order” can be computed in \(O(N)\) time using the principle of merge sort (you can use C++’s std::merge or Python’s heapq.merge), so the original problem can be solved in \(O(N (\log N + \log M))\) time. (If you are using a fast language, you could afford to spend \(O(N \log N)\) time to sort the \(2N\) values every time, making the overall time \(O(N \log N \log M)\), which should just about be feasible.)

Sample Implementation (Python):

import heapq

N, K = map(int, input().split())
X, Y = [], []
for _ in range(N):
    x, y = map(int, input().split())
    X.append(x)
    Y.append(y)
X.sort()
Y.sort()

def cost(z, l):
    p = list(heapq.merge(z, map(lambda x: x - l, z)))[N]
    c = 0
    for x in z:
        if x < p:
            c += p - x
        elif x > p + l:
            c += x - (p + l)
    return c

lo, hi = -1, 10**9
while hi - lo > 1:
    mi = (lo + hi) // 2
    if cost(X, mi) + cost(Y, mi) <= K:
        hi = mi
    else:
        lo = mi
print(hi)