(Analysis by Mark Chen)

For any set of points P, the minimum area of an enclosing rectangle with sides parallel to the x and y axes is equal to the product of the difference between the largest and smallest x-coordinates in P and the difference between the largest and smallest y-coordinates in P. For a proof sketch, note that the projection of any such rectangle onto the x-axis must contain the interval [min x-coordinate in P, max x-coordinate in P], or it cannot have contained all points in P.

Now consider the four largest x-coordinates (allowing repeats) over all points in P. After we remove at most 3 points, the resulting largest x-coordinate must be one of these four values. Similarly, after we remove at most 3 points, there are only four possible values for the smallest x-coordinate, the largest y-coordinate, and the smallest y-coordinate.

Since there are just four candidates for each side of the new rectangle, there are at most $4^4 = 256$ possible rectangles that could result from removing 3 points! For each candidate rectangle, we iterate through all points in P to count how many points lie outside of it. If this count is less than or equal to 3, we have a valid rectangle, and should compute its area. The final answer is the minimum of all valid rectangle areas.

Here is Travis Hance's code:

#include <cstdio>
#include <vector>
#include <algorithm>

#define NMAX 100000

int n;
long long x[NMAX];
long long y[NMAX];
#define infinite 1000000000

struct Analysis {
long long area;
std::vector<std::vector<int> > borders;
};

Analysis analyze(std::vector<int> indicesToSkip) {
long long minX = infinite, minY = infinite, maxX = -infinite, maxY = -infinite;
for (int i = 0; i < n; i++) {
bool skip = false;
for (int j = 0; j < indicesToSkip.size(); j++) {
if (indicesToSkip[j] == i) {
skip = true;
}
}

if (skip) continue;

minX = std::min(minX, x[i]);
maxX = std::max(maxX, x[i]);
minY = std::min(minY, y[i]);
maxY = std::max(maxY, y[i]);
}

Analysis a;
a.area = (maxX - minX) * (maxY - minY);

std::vector<int> up, down, left, right;

for (int i = 0; i < n; i++) {
bool skip = false;
for (int j = 0; j < indicesToSkip.size(); j++) {
if (indicesToSkip[j] == i) {
skip = true;
}
}

if (skip) continue;

if (x[i] == minX) left.push_back(i);
if (x[i] == maxX) right.push_back(i);
if (y[i] == minY) up.push_back(i);
if (y[i] == maxY) down.push_back(i);
}

if (up.size() <= 3) a.borders.push_back(up);
if (down.size() <= 3) a.borders.push_back(down);
if (left.size() <= 3) a.borders.push_back(left);
if (right.size() <= 3) a.borders.push_back(right);

return a;
}

int main() {
freopen("reduce.in", "r", stdin);
freopen("reduce.out", "w", stdout);

scanf("%d", &n);
for (int i = 0; i < n; i++) {
scanf("%lld", &x[i]);
scanf("%lld", &y[i]);
}

Analysis a = analyze(std::vector<int>());
long long bestArea = a.area;

for (std::vector<int> pointsOnBorder : a.borders) {
Analysis smallerAnalysis = analyze(pointsOnBorder);
bestArea = std::min(bestArea, smallerAnalysis.area);
for (std::vector<int> pointsOnBorder2 : smallerAnalysis.borders) {
if (pointsOnBorder2.size() + pointsOnBorder.size() <= 3) {
for (int p : pointsOnBorder) {
pointsOnBorder2.push_back(p);
}
Analysis analysis3 = analyze(pointsOnBorder2);
bestArea = std::min(bestArea, analysis3.area);
for (std::vector<int> pointsOnBorder3 : analysis3.borders) {
if (pointsOnBorder2.size() + pointsOnBorder3.size() <= 3) {
for (int p : pointsOnBorder2) {
pointsOnBorder3.push_back(p);
}
Analysis analysis4 = analyze(pointsOnBorder3);
bestArea = std::min(bestArea, analysis4.area);
}
}
}
}
}

printf("%lld\n", bestArea);
}