#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MOD = 1E9 + 7;
const int INF = 1E9; const ll INFLL = 1E18;
const int MAX = 2001;
const int LOG = 11;
int lg[MAX];
struct SparseTable {
int n;
vector<vector<int>> sparse;
SparseTable(int n) {
this->n = n;
sparse.assign(n, vector<int>(LOG));
}
void build(vector<int>& a) {
for(int i = 0; i < n; i++) {
sparse[i][0] = a[i];
}
for(int i = 1; i < LOG; i++) {
for(int j = 0; j + (1 << i) - 1 < n; j++) {
sparse[j][i] = max(sparse[j][i - 1], sparse[j + (1 << (i - 1))][i - 1]);
}
}
}
int query(int l, int r) {
int k = lg[r - l + 1];
return max(sparse[l][k], sparse[r - (1 << k) + 1][k]);
}
};
int biggest_stadium(int N, vector<vector<int>> F) {
lg[1] = 0;
for(int i = 2; i < MAX; i++) {
lg[i] = lg[i / 2] + 1;
}
vector<vector<int>> treeUp(N, vector<int>(N, -1));
vector<vector<int>> treeDown(N, vector<int>(N, N));
vector<SparseTable> sUp(N, SparseTable(N));
vector<SparseTable> sDown(N, SparseTable(N));
vector<set<int>> treeRows(N);
for(int i = 0; i < N; i++) {
for(int j = 0; j < N; j++) {
if(i) {
treeUp[i][j] = treeUp[i - 1][j];
}
if(F[i][j]) {
treeUp[i][j] = i;
treeRows[i].insert(j);
}
}
}
for(int i = N - 1; i >= 0; i--) {
for(int j = 0; j < N; j++) {
if(i + 1 < N) {
treeDown[i][j] = treeDown[i + 1][j];
}
if(F[i][j]) {
treeDown[i][j] = i;
}
}
}
for(int i = 0; i < N; i++) {
for(int j = 0; j < N; j++) {
treeDown[i][j] *= -1;
}
sUp[i].build(treeUp[i]);
sDown[i].build(treeDown[i]);
}
vector<map<array<int, 4>, int>> width(N);
array<int, 4> v = {0, 0, 0, 0};
for(int i = 0; i < N; i++) {
for(int j = 0; j < N; j++) {
int k = j;
while(k < N && !F[i][k]) {
k++;
}
if(!F[i][j]) {
int wid = k - j - 1;
int x1; int y1; int x2; int y2;
y1 = j; y2 = k - 1;
x1 = sUp[i].query(y1, y2) + 1;
x2 = -sDown[i].query(y1, y2) - 1;
width[wid][{x1, y1, x2, y2}] = (y2 - y1 + 1) * (x2 - x1 + 1);
}
j = k;
}
}
int ans = 0;
for(int i = N - 1; i >= 0; i--) {
for(pair<array<int, 4>, int> j : width[i]) {
ans = max(ans, j.second);
int x1 = j.first[0]; int y1 = j.first[1]; int x2 = j.first[2]; int y2 = j.first[3];
if(x1) {
vector<int> pts(1, y1 - 1);
set<int>::iterator it = treeRows[x1 - 1].lower_bound(y1);
while(it != treeRows[x1 - 1].end() && *it <= y2) {
pts.push_back(*it);
++it;
}
pts.push_back(y2 + 1);
for(int k = 0; k + 1 < pts.size(); k++) {
if(pts[k + 1] > pts[k] + 1) {
int nx1; int ny1; int nx2; int ny2;
ny1 = pts[k] + 1; ny2 = pts[k + 1] - 1;
nx1 = sUp[x1 - 1].query(ny1, ny2) + 1;
if(x2 + 1 < N) {
nx2 = -sDown[x2 + 1].query(ny1, ny2) - 1;
} else {
nx2 = x2;
}
int wid = ny2 - ny1;
int addArea = (ny2 - ny1 + 1) * (nx2 - nx1 + 1) - (ny2 - ny1 + 1) * (x2 - x1 + 1);
width[wid][{nx1, ny1, nx2, ny2}] = max(width[wid][{nx1, ny1, nx2, ny2}], j.second + addArea);
}
}
}
if(x2 + 1 < N) {
vector<int> pts(1, y1 - 1);
set<int>::iterator it = treeRows[x2 + 1].lower_bound(y1);
while(it != treeRows[x2 + 1].end() && *it <= y2) {
pts.push_back(*it);
++it;
}
pts.push_back(y2 + 1);
for(int k = 0; k + 1 < pts.size(); k++) {
if(pts[k + 1] > pts[k] + 1) {
int nx1; int ny1; int nx2; int ny2;
ny1 = pts[k] + 1; ny2 = pts[k + 1] - 1;
if(x1) {
nx1 = sUp[x1 - 1].query(ny1, ny2) + 1;
} else {
nx1 = x1;
}
nx2 = -sDown[x2 + 1].query(ny1, ny2) - 1;
int wid = ny2 - ny1;
int addArea = (ny2 - ny1 + 1) * (nx2 - nx1 + 1) - (ny2 - ny1 + 1) * (x2 - x1 + 1);
width[wid][{nx1, ny1, nx2, ny2}] = max(width[wid][{nx1, ny1, nx2, ny2}], j.second + addArea);
}
}
}
}
}
return ans;
}
