#include <bits/stdc++.h>
using namespace std;
template <typename T, class F = function<T(const T&, const T&)>>
class SparseTable {
public:
int n;
vector<vector<T>> mat;
F func;
SparseTable(){};
SparseTable(const vector<T>& a, const F& f) : func(f) {
n = static_cast<int>(a.size());
int max_log = 32 - __builtin_clz(n);
mat.resize(max_log);
mat[0] = a;
for (int j = 1; j < max_log; j++) {
mat[j].resize(n - (1 << j) + 1);
for (int i = 0; i <= n - (1 << j); i++) {
mat[j][i] = func(mat[j - 1][i], mat[j - 1][i + (1 << (j - 1))]);
}
}
}
T get(int from, int to) const {
assert(0 <= from && from <= to && to <= n - 1);
int lg = 32 - __builtin_clz(to - from + 1) - 1;
return func(mat[lg][from], mat[lg][to - (1 << lg) + 1]);
}
};
int biggest_stadium(int n, vector<vector<int>> F) {
vector<vector<int>> row(n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
if (F[i][j]) {
row[i].push_back(j);
}
}
}
vector<SparseTable<int>> Lrec(n);
vector<SparseTable<int>> Hrec(n);
vector<int> cur(n, -1);
for (int r = 0; r < n; r++) {
for (int c = 0; c < n; c++) {
if (F[r][c]) cur[c] = r;
}
Hrec[r] = SparseTable<int>(cur, [&](int i, int j){ return max(i, j); });
}
cur = vector<int>(n, n);
for (int r = n - 1; r >= 0; r--) {
for (int c = 0; c < n; c++) {
if (F[r][c]) cur[c] = r;
}
Lrec[r] = SparseTable<int>(cur, [&](int i, int j){ return min(i, j); });
}
auto get_rec = [&](int y, int Lx, int Rx) {
int Ly = Hrec[y].get(Lx, Rx) + 1;
int Ry = Lrec[y].get(Lx, Rx) - 1;
return make_tuple(Lx, Rx, Ly, Ry);
};
int ans = 0;
map<tuple<int, int, int, int>, int> dp;
function<int(int, int, int, int)> dfs = [&](int Lx, int Rx, int Ly, int Ry) {
if (dp.find(make_tuple(Lx, Rx, Ly, Ry)) != dp.end()) return dp[make_tuple(Lx, Rx, Ly, Ry)];
int ret = (Rx - Lx + 1) * (Ry - Ly + 1);
auto solve_row = [&](int y) {
vector<int> que;
for (auto it = lower_bound(row[y].begin(), row[y].end(), Lx); it != row[y].end() && *it <= Rx; it++) {
que.push_back(*it);
}
assert(que.size());
if (que[0] > Lx) {
auto [lx, rx, ly, ry] = get_rec(y, Lx, que[0] - 1);
int add = (Ry - Ly + 1) * (Rx - Lx - (rx - lx));
ret = max(ret, add + dfs(lx, rx, ly, ry));
}
for (int i = 0; i < (int) que.size() - 1; i++) {
if (que[i + 1] - que[i] == 1) continue;
auto [lx, rx, ly, ry] = get_rec(y, que[i] + 1, que[i + 1] - 1);
int add = (Ry - Ly + 1) * (Rx - Lx - (rx - lx));
ret = max(ret, add + dfs(lx, rx, ly, ry));
}
if (que.back() < Rx) {
auto [lx, rx, ly, ry] = get_rec(y, que.back() + 1, Rx);
int add = (Ry - Ly + 1) * (Rx - Lx - (rx - lx));
ret = max(ret, add + dfs(lx, rx, ly, ry));
}
};
if (Ly > 0) {
solve_row(Ly - 1);
}
if (Ry < n - 1) {
solve_row(Ry + 1);
}
return dp[make_tuple(Lx, Rx, Ly, Ry)] = ret;
};
for (int r = 0; r < n; r++) {
vector<int> cur = {-1};
for (int c : row[r]) cur.push_back(c);
cur.push_back(n);
int sz = cur.size();
for (int i = 0; i < sz - 1; i++) {
int L = cur[i], R = cur[i + 1];
if (R == L + 1) continue;
auto [Lx, Rx, Ly, Ry] = get_rec(r, L + 1, R - 1);
assert(Lx == L + 1 && Rx == R - 1 && Ly <= r && Ry >= r);
ans = max(ans, dfs(Lx, Rx, Ly, Ry));
}
}
return ans;
}
