We define a function \(f(A, B, C, D)\) as follows. (Notice the presence/absence of equal sign in each inequality sign.)

• \(f(A, B, C, D) =\) (the number of black squares among those squares \((h, w)\) with \(A \leq h \lt C\) and \(B \leq w \lt D\)).

The sought value is \(f(A, B, C + 1, D + 1)\). (Note the shift by \(1\).)
While one can directly implement an algorithm that directly evaluates \(f\), the implementation requires a bunch of caseworks because there are four arguments, so we try to avoid it.

What we want is the (d) part:

By exploiting a technique on cumulative sums, we can express it as follows:

\[ \begin{aligned} &f(A, B, C, D) \\ &=(\text{The number of black squares in }(d)) \\ &= (\text{The number of black squares in }(a), (b), (c), (d)) \\ &- (\text{The number of black squares in }(a), (b)) \\ &- (\text{The number of black squares in }(a), (c)) \\ &+ (\text{The number of black squares in }(a)) \\ &= f(0, 0, C, D) - f(0, 0, C, B) - f(0, 0, A, D) + f(0, 0, A, B) \end{aligned} \]

With this deformation, we can constrain the arguments of \(f\) so that \(A=B=0\). That is, if we define

• \(g(H, W) =\) (the number of black squares among those squares \((h, w)\) with \(0 \leq h \lt H\) and \(0 \leq w \lt W\)),

it holds that

\[ f(A,B,C,D) = g(C,D) - g(C, B) - g(A, D) + g(A, B), \]

so all that left is to implement \(g\), which has only two arguments.

Now we consider how to evaluate \(g(H, W)\). A naive implementation would cost too much time, so we want to make use of the periodicity.
For example, when \(H = 8, W = 7\), and \(N = 3\), one can divide the grid by period as follows:

Since the grid has \((3 \times 3)\) periodic pattern, the sub-rectangles are classified into the following four types, all of which coincide with or are smaller than \(N \times N\):

• \(4\) rectangles of type A whose pattern coincides with the top-left \(3 \times 3\) region
• \(2\) rectangles of type B whose pattern coincides with the top-left \(2 \times 3\) region
• \(2\) rectangles of type C whose pattern coincides with the top-left \(3 \times 1\) region
• \(1\) rectangles of type D whose pattern coincides with the top-left \(2 \times 1\) region

This holds for a general \(H, W, N\), so the problem can be simplified if we know the number of rectangles of the four types and their dimensions. These can be obtained as the quotients and remainders when \(H\) and \(W\) are divided by \(N\).

Therefore, by precomputing \(g(i, j)\) for all \(0 \leq i \leq N, 0 \leq j \leq N\), the problem can be solved in \(\mathrm{O}(1)\) time per query.
The percomputation can be done by applying the technique used for the first figure again, in order to obtain the following recurrence relation:

\[ \begin{aligned} g(i, j) &= (1\text{ if } ((i-1,j-1)\text{ is a black square}) \text{ else } 0) \\ &+ g(i - 1, j) + g(i, j - 1) - g(i - 1, j - 1). \end{aligned} \]

Thus, the precalculation can be done in a total of \(\mathrm{O}(N^2)\) time. (If \(i = 0\) or \(j = 0\), let \(g(i, j) = 0\).)
Therefore, the problem has been solved in a total of \(\mathrm{O}(N^2 + Q)\) time, which is fast enough.

• Sample code (C++)

#include <iostream>
#include <string>
#include <vector>
using namespace std;

int N, Q, precalc[1010][1010];

long long g(int H, int W) {
  if (H <= N and W <= N) return precalc[H][W];
  int Hq = H / N, Hr = H % N;
  int Wq = W / N, Wr = W % N;
  long long ret = 0;
  ret += g(N, N) * Hq * Wq;
  ret += g(Hr, N) * Wq;
  ret += g(N, Wr) * Hq;
  ret += g(Hr, Wr);
  return ret;
}

long long f(int A, int B, int C, int D) {
  return g(C, D) - g(A, D) - g(C, B) + g(A, B);
}

int main() {
  cin >> N >> Q;
  vector<string> S(N);
  for (auto& s : S) cin >> s;
  for (int i = 1; i <= N; i++) {
    for (int j = 1; j <= N; j++) {
      precalc[i][j] += S[i - 1][j - 1] == 'B';
      precalc[i][j] += precalc[i - 1][j];
      precalc[i][j] += precalc[i][j - 1];
      precalc[i][j] -= precalc[i - 1][j - 1];
    }
  }
  while (Q--) {
    int A, B, C, D;
    cin >> A >> B >> C >> D;
    cout << f(A, B, C + 1, D + 1) << "\n";
  }
}