Problem: Problem F and F2 are the same problem, but with different constraints and time limits.

We have a board divided into N horizontal rows and N vertical columns of square cells. The cell at the i-th row from the top and the j-th column from the left is called Cell (i,j). Each cell is either empty or occupied by an obstacle. Also, each empty cell has a digit written on it. If A_{i,j}= `1`, `2`, ..., or `9`, Cell (i,j) is empty and the digit A_{i,j} is written on it. If A_{i,j}= `#`, Cell (i,j) is occupied by an obstacle.

Cell Y is reachable from cell X when the following conditions are all met:

* Cells X and Y are different.
* Cells X and Y are both empty.
* One can reach from Cell X to Cell Y by repeatedly moving right or down to an adjacent empty cell.



Consider all pairs of cells (X,Y) such that cell Y is reachable from cell X. Find the sum of the products of the digits written on cell X and cell Y for all of those pairs.

Constraints

* 1 \leq N \leq 500
* A_{i,j} is one of the following characters: `1`, `2`, ... `9` and `#`.

Input

Input is given from Standard Input in the following format:


N
A_{1,1}A_{1,2}...A_{1,N}
A_{2,1}A_{2,2}...A_{2,N}
:
A_{N,1}A_{N,2}...A_{N,N}


Output

Print the sum of the products of the digits written on cell X and cell Y for all pairs (X,Y) such that cell Y is reachable from cell X.

Examples

Input

2
11
11


Output

5


Input

4
1111
11#1
1#11
1111


Output

47


Input

10
76##63##3#
8445669721
75#9542133
3#285##445
749632##89
2458##9515
5952578#77
1#3#44196#
4355#99#1#
298#63587


Output

36065


Input

10
4177143673
7#########
5#1716155#
6#4#####5#
2#3#597#6#
6#9#8#3#5#
5#2#899#9#
1#6#####6#
6#5359657#
5#########


Output

6525

Solution (in CPP):

Solution:
// includes {{{
#include<iostream>
#include<iomanip>
#include<algorithm>
#include<vector>
#include<stack>
#include<queue>
#include<deque>
#include<map>
#include<set>
#include<tuple>
#include<cmath>
#include<random>
#include<cassert>
// }}}
// #include<bits/stdc++.h>
using namespace std;
using ll = long long;

/*

   checklist
   - min vs max
   - ymin + 1
   - for range

 */

// O(N^2 log N)
// 想定解法と同じ

template<class T, class U> inline void smin(T &a, const U &b) { a = a < (T) b ? a : b; }
template<class T, class U> inline void smax(T &a, const U &b) { a = a < (T) b ? b : a; }

const int N = 1501;
string s[N]; // y, x
int maxX[N][N], maxY[N][N];
int minX[N][N], minY[N][N];
const int inf = 1e9;
// 9 * N * N < inf
// ans max = O(9 * 9 * N^4)

int get(int y, int x, int f) {
  if(f) swap(y, x);
  return s[y][x];
}

int getMinX(int y, int x, int f) {
  if(f) swap(y, x);
  return f ? minY[y][x] : minX[y][x];
}
int getMinY(int y, int x, int f) {
  if(f) swap(y, x);
  return f ? minX[y][x] : minY[y][x];
}
int getMaxX(int y, int x, int f) {
  if(f) swap(y, x);
  return f ? maxY[y][x] : maxX[y][x];
}
int getMaxY(int y, int x, int f) {
  if(f) swap(y, x);
  return f ? maxX[y][x] : maxY[y][x];
}

template<class T, class U = ll>
U range2D(T &val, int y1, int x1, int y2, int x2, int h, int w) {
  if(y1 < 0) y1 = 0;
  if(x1 < 0) x1 = 0;
  if(y2 >= h) y2 = h - 1;
  if(x2 >= w) x2 = w - 1;
  if(y1 > y2 || x1 > x2) return U(0);
  U res = val[y2][x2];
  if(y1-1 >= 0) res -= val[y1-1][x2];
  if(x1-1 >= 0) res -= val[y2][x1-1];
  if(y1-1 >= 0 && x1-1 >= 0) res += val[y1-1][x1-1];
  return res;
}

template<class T>
void accum2D(T &val, int h, int w) {
  for(int i = 1; i < h; i++) for(int j = 0; j < w; j++) val[i][j] += val[i-1][j];
  for(int i = 0; i < h; i++) for(int j = 1; j < w; j++) val[i][j] += val[i][j-1];
}

// DAQ
ll solve(int y1, int x1, int y2, int x2, int f = 0) {
  int h = y2 - y1 + 1;
  int w = x2 - x1 + 1;

  // 自明 {{{
  if(h == 1 || w == 1) {
    if(h == 1) {
      int noe = 0;
      ll res = 0;
      for(int x = x1; x <= x2; x++) {
        int eca = get(y1, x, f);
        if(eca) {
          res += (ll) eca * noe;
          noe += eca;
        } else noe = 0;
      }
      return res;
    } else {
      return solve(x1, y1, x2, y2, !f);
    }
  }
  // }}}

  // 上下に分けたものについてのみ解く
  // w <= h

  int ymid = (y1 + y2) >> 1;
  // ymid自体は上側

  ll res = 0;
  // 上のほうが大きい
  // 渡すときにswapする
  res += solve(x1, y1, x2, ymid, !f);
  res += solve(x1, ymid + 1, x2, y2, !f);
  int h0 = ymid - y1 + 1;
  int h1 = h - h0;

  // ローカルな関数 {{{

  vector<vector<int>> minMeet(w, vector<int>(w, inf)); // y座標の最小値
  vector<vector<int>> L(h1, vector<int>(w, inf));
  vector<vector<int>> R(h1, vector<int>(w, -inf));

  vector<vector<int>> sumByYL(h1, vector<int>(w, 0));
  vector<vector<int>> sumByYR(h1, vector<int>(w, 0));
  vector<vector<int>> sumByLR(w, vector<int>(w, 0));

  vector<int> validSum(h1);

  auto getMinMax = [&](int a, int b) {
    return min(getMaxY(ymid + 1, x1 + a, f), getMaxY(ymid + 1, x1 + b, f));
  };

  auto getBothReachable = [&](int a, int b, int l) {
    smin(l, getMinMax(a, b) - (ymid + 1)); // これは必要で，サンプルケース4がhack
    if(minMeet[a][b] > l) return 0;
    int res = validSum[l];
    res -= range2D(sumByYL, 0, a + 1, l, w - 1, h1, w);
    res -= range2D(sumByYR, 0, 0, l, b - 1, h1, w);
    res += range2D(sumByLR, a + 1, 0, w - 1, b - 1, w, w);
    return res;
  };

  // }}}
  
  // 前計算1 {{{

  for(int j = 0; j < w; j++) {
    if(get(ymid + 1, x1 + j, f)) {
      L[0][j] = j;
      R[0][j] = j;
    }
  }
  for(int i = 0; i < h1; i++) {
    for(int j = 0; j < w; j++) {
      if(get(ymid + 1 + i, x1 + j, f)) {
        if(i) smin(L[i][j], L[i-1][j]);
        if(j) smin(L[i][j], L[i][j-1]);
        if(i) smax(R[i][j], R[i-1][j]);
        if(j) smax(R[i][j], R[i][j-1]);
        // 値が決定
        if(L[i][j] <= R[i][j]) {
          smin(minMeet[L[i][j]][R[i][j]], i);
        }
      }
    }
  }

  // 累積min
  for(int i = 0; i < w; i++) for(int j = w - 2; j >= 0; j--) smin(minMeet[i][j], minMeet[i][j+1]);
  for(int i = 1; i < w; i++) for(int j = w - 1; j >= 0; j--) smin(minMeet[i][j], minMeet[i-1][j]);

  for(int i = 0; i < w; i++) for(int j = w - 1; j >= 0; j--) {
    if(getMinMax(i, j) < minMeet[i][j]) {
      minMeet[i][j] = inf;
    }
  }
  
  // }}}

  // 前計算2 {{{

  for(int i = 0; i < h1; i++) for(int j = 0; j < w; j++) {
    if(L[i][j] <= R[i][j]) {
      int ecas = get(ymid + 1 + i, x1 + j, f); // != 0
      sumByYL[i][L[i][j]] += ecas;
      sumByYR[i][R[i][j]] += ecas;
      sumByLR[L[i][j]][R[i][j]] += ecas;
      validSum[i] += ecas;
    }
  }

  for(int i = 1; i < h1; i++) validSum[i] += validSum[i-1];

  accum2D(sumByYL, h1, w);
  accum2D(sumByYR, h1, w);
  accum2D(sumByLR, w, w);

  // }}}

  vector<vector<int>> UL(h0, vector<int>(w, inf));
  vector<vector<int>> UR(h0, vector<int>(w, -inf));

  // 前計算3 {{{

  for(int j = 0; j < w; j++) {
    if(get(ymid, x1 + j, f)) {
      UL[h0 - 1][j] = j;
      UR[h0 - 1][j] = j;
    }
  }

  for(int i = h0 - 1; i >= 0; i--) for(int j = w - 1; j >= 0; j--) {
    if(get(y1 + i, x1 + j, f)) {
      if(i+1 < h0) smin(UL[i][j], UL[i+1][j]);
      if(j+1 < w) smin(UL[i][j], UL[i][j+1]);
      if(i+1 < h0) smax(UR[i][j], UR[i+1][j]);
      if(j+1 < w) smax(UR[i][j], UR[i][j+1]);
    }
  }

  // }}}


  vector<int> qina(w); // i+1からは取れないが，iからは取れるものの sum
  for(int i = 0; i < h0; i++) {
    int l = 0, r = 0; // [l, r)
    int sum = 0;
    stack<pair<int, int>> tapris; // index, height
    for(int j = 0; j < w; j++) if(UL[i][j] <= UR[i][j]) {
      // U(y1+i, x1+j) は必ず D(0, _) のどこかにたどり着く
      // add
      for(;r < w && r <= UR[i][j]; r++) {
        if(getMinY(ymid, x1 + r, f) > y1 + i) {
          // 無効値．D(0, x1 + r) は U(y1+i, _) にたどり着かない
          qina[r] = 0;
          continue;
        }
        qina[r] = getBothReachable(r, r, h1 - 1);
        assert(qina[r] >= 0);
        sum += qina[r];
        int hei = min(getMaxY(ymid + 1, x1 + r, f) - (ymid + 1), h1 - 1);
        int last = -1;
        auto f = [&]() {
          int idx, k;
          tie(idx, k) = tapris.top();
          if(qina[idx] == -1) return; // 削除済みの値
          sum -= qina[idx];
          qina[idx] -= getBothReachable(idx, r, h1 - 1) - (last == -1 ? 0 : getBothReachable(idx, r, last));
          // assert(qina[idx] >= 0);
          sum += qina[idx];
          last = k;
        };
        while(tapris.size() && tapris.top().second <= hei) {
          f();
          tapris.pop();
        }
        if(tapris.size()) {
          f();
        }
        tapris.emplace(r, hei);
      }
      // rem
      while(l < w && l < UL[i][j]) {
        sum -= qina[l];
        qina[l] = -1;
        l++;
      }
      //
      res += (ll) sum * get(y1 + i, x1 + j, f);
    }
  }

  // if(f) swap(y1, x1), swap(y2, x2);
  return res;
}

int main() {
  std::ios::sync_with_stdio(false), std::cin.tie(0);
  int n; cin >> n;
  for(int i = 0; i < n; i++) {
    cin >> s[i];
    for(char &c : s[i]) c = c == '#' ? 0 : c - '0';
  }

  // グローバルな前計算 {{{
  // maxX
  for(int i = n - 1; i >= 0; i--) for(int j = n - 1; j >= 0; j--) {
    if(s[i][j]) {
      maxX[i][j] = j;
      if(i + 1 < n) smax(maxX[i][j], maxX[i + 1][j]);
      if(j + 1 < n) smax(maxX[i][j], maxX[i][j + 1]);
    } else maxX[i][j] = -inf;
  }
  // maxY
  for(int i = n - 1; i >= 0; i--) for(int j = n - 1; j >= 0; j--) {
    if(s[j][i]) {
      maxY[j][i] = j;
      if(i + 1 < n) smax(maxY[j][i], maxY[j][i + 1]);
      if(j + 1 < n) smax(maxY[j][i], maxY[j + 1][i]);
    } else maxY[j][i] = -inf;
  }
  // minX
  for(int i = 1; i < n; i++) for(int j = 1; j < n; j++) {
    if(s[i][j]) {
      minX[i][j] = j;
      if(i - 1 >= 0) smin(minX[i][j], minX[i - 1][j]);
      if(j - 1 >= 0) smin(minX[i][j], minX[i][j - 1]);
    } else minX[i][j] = inf;
  }
  // minY
  for(int i = 1; i < n; i++) for(int j = 1; j < n; j++) {
    if(s[j][i]) {
      minY[j][i] = j;
      if(i - 1 >= 0) smin(minY[j][i], minY[j][i - 1]);
      if(j - 1 >= 0) smin(minY[j][i], minY[j - 1][i]);
    } else minY[j][i] = inf;
  }
  // }}}

  cout << solve(0, 0, n - 1, n - 1) << "\n";
  return 0;
}