関数 \(f(A, B, C, D)\) を次のように定義します。(不等号のイコールの有無に注意してください)

• \(f(A, B, C, D) =\) (\(A \leq h \lt C, B \leq w \lt D\) を満たす範囲にあるマス \((h, w)\) のうち黒マスであるマスの個数)

求めたい値は \(f(A, B, C + 1, D + 1)\) です (\(1\) ズレるのに注意してください)
この \(f\) を計算するアルゴリズムをそのまま実装しても解くことが出来ますが、このままでは引数が \(4\) 個あるせいで細かい場合分けが増えて面倒なので、場合分けを減らすことを試みます。

求めたいものは下の図の (d) の部分です。

ここで累積和に関するテクニックを利用すると、次のような式変形が可能です。

\[ \begin{aligned} &f(A, B, C, D) \\ &=((d) に含まれる黒マスの個数) \\ &= ((a), (b), (c), (d) に含まれる黒マスの個数) \\ &- ((a), (b) に含まれる黒マスの個数) \\ &- ((a), (c) に含まれる黒マスの個数) \\ &+ ((a) に含まれる黒マスの個数) \\ &= f(0, 0, C, D) - f(0, 0, C, B) - f(0, 0, A, D) + f(0, 0, A, B) \end{aligned} \]

このように変形すると、\(f\) は \(A=B=0\) の場合のみを考えればよいことになります。つまり、

• \(g(H, W) =\) (\(0 \leq h \lt H, 0 \leq w \lt W\) を満たす範囲にあるマス \((h, w)\) のうち黒マスであるマスの個数)

と置いたとき

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

が成り立つので、\(g\) を計算する関数を書けばよく、引数が \(2\) 個に減っています。

さて、\(g(H, W)\) を計算する方法を考えます。これは愚直にやると大変な時間がかかってしまうので、うまく周期性を利用することを目指ししょう。
例えば \(H = 8, W = 7, N = 3\) の場合、グリッドを周期で切ると下図のようになります。

グリッドは \(3 \times 3\) 周期で同じ模様が現れるので、上の図に含まれる分割された長方形領域は次の \(4\) 種類のみで、これらは全て \(N \times N\) と一致するかそれよりも小さいです。

• 左上 \(3 \times 3\) と一致する模様の長方形 A が \(4\) 個
• 左上 \(2 \times 3\) と一致する模様の長方形 B が \(2\) 個
• 左上 \(3 \times 1\) と一致する模様の長方形 C が \(2\) 個
• 左上 \(2 \times 1\) と一致する模様の長方形 D が \(1\) 個

\(H, W, N\) が一般の場合でもこのような性質は成り立つため、4 種類の長方形の個数と種類が分かれば問題を単純化することが出来ます。これは \(H, W\) を \(N\) で割った商と余りを利用することで個数と種類を計算することが出来ます。

以上より、\(g(i, j)\) をあらかじめ \(0 \leq i \leq N, 0 \leq j \leq N\) の範囲で前計算しておけば、クエリを \(\mathrm{O}(1)\) で処理することが出来るようになりました。
\(g(i, j)\) の前計算は \(1\) 枚目の画像で用いたテクニックをここでも利用することで

\[ \begin{aligned} g(i, j) &= (1\text{ if } ((i-1,j-1) が黒マス) \text{ else } 0) \\ &+ g(i - 1, j) + g(i, j - 1) - g(i - 1, j - 1) \end{aligned} \]

という漸化式を導出することが出来るので、ここから \(\mathrm{O}(N^2)\) で前計算することが出来ます。 (\(i = 0\) または \(j = 0\) の場合は \(g(i, j) = 0\) とする)
以上よりこの問題を \(\mathrm{O}(N^2 + Q)\) で解くことが出来ました。これは十分高速です。

• 実装例(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";
  }
}