Problem:
 Score : 700 points    Problem Statement Nuske has a grid with N rows and M columns of squares. The rows are numbered 1 through N from top to bottom, and the columns are numbered 1 through M from left to right.
Each square in the grid is painted in either blue or white. If S_{i,j} is 1 , the square at the i -th row and j -th column is blue; if S_{i,j} is 0 , the square is white.
For every pair of two blue square a and b , there is at most one path that starts from a , repeatedly proceeds to an adjacent (side by side) blue square and finally reaches b , without traversing the same square more than once.  Phantom Thnook, Nuske's eternal rival, gives Q queries to Nuske. The i -th query consists of four integers x_{i,1} , y_{i,1} , x_{i,2} and y_{i,2} and asks him the following: when the rectangular region of the grid bounded by (and including) the x_{i,1} -th row, x_{i,2} -th row, y_{i,1} -th column and y_{i,2} -th column is cut out, how many connected components consisting of blue squares there are in the region?  Process all the queries.      Constraints  1 ≤ N,M ≤ 2000  1 ≤ Q ≤ 200000  S_{i,j} is either 0 or 1 .  S_{i,j} satisfies the condition explained in the statement.  1 ≤ x_{i,1} ≤ x_{i,2} ≤ N(1 ≤ i ≤ Q)  1 ≤ y_{i,1} ≤ y_{i,2} ≤ M(1 ≤ i ≤ Q)         Input The input is given from Standard Input in the following format:  N  M  Q  S_{1,1} .. S_{1,M} : S_{N,1} .. S_{N,M}  x_{1,1}  y_{i,1}  x_{i,2}  y_{i,2} : x_{Q,1}  y_{Q,1}  x_{Q,2}  y_{Q,2}       Output For each query, print the number of the connected components consisting of blue squares in the region.        Sample Input 1 3 4 4
1101
0110
1101
1 1 3 4
1 1 3 1
2 2 3 4
1 2 2 4      Sample Output 1 3
2
2
2   In the first query, the whole grid is specified. There are three components consisting of blue squares, and thus 3 should be printed.  In the second query, the region within the red frame is specified. There are two components consisting of blue squares, and thus 2 should be printed.
Note that squares that belong to the same component in the original grid may belong to different components.       Sample Input 2 5 5 6
11010
01110
10101
11101
01010
1 1 5 5
1 2 4 5
2 3 3 4
3 3 3 3
3 1 3 5
1 1 3 4      Sample Output 2 3
2
1
1
3
2

Solution (Python):
import sys
input = sys.stdin.readline

N, M, Q = map(int, input().split())
S = [input().rstrip() for _ in range(N)]
cum = [[0] * (M + 1) for _ in range(N + 1)]
for i in range(N):
    for j in range(M):
        cum[i+1][j+1] = int(S[i][j]) + cum[i][j+1] + cum[i+1][j] - cum[i][j]
cumh = [[0] * M for _ in range(N + 1)]
for i in range(N):
    for j in range(M-1):
        cumh[i+1][j+1] = cumh[i][j+1] + cumh[i+1][j] - cumh[i][j]
        if S[i][j] == S[i][j+1] == '1':
            cumh[i+1][j+1] += 1
cumv = [[0] * (M + 1) for _ in range(N)]
for i in range(N-1):
    for j in range(M):
        cumv[i+1][j+1] = cumv[i][j+1] + cumv[i+1][j] - cumv[i][j]
        if S[i][j] == S[i+1][j] == '1':
            cumv[i+1][j+1] += 1
for _ in range(Q):
    x1, y1, x2, y2 = map(lambda x: int(x) - 1, input().split())
    a = cum[x2+1][y2+1] - cum[x2+1][y1] - cum[x1][y2+1] + cum[x1][y1]
    b = cumh[x2+1][y2] - cumh[x1][y2] - cumh[x2+1][y1] + cumh[x1][y1]
    c = cumv[x2][y2+1] - cumv[x2][y1] - cumv[x1][y2+1] + cumv[x1][y1]
    print(a - b - c)

Solution (Python):
import numpy as np

N, M, Q = map(int, input().split())
field = np.array([list(map(int, input())) for _ in range(N)], dtype = np.int64)
xyxys = [tuple(map(int, input().split())) for _ in range(Q)]

above = np.zeros((N+1, M+1), dtype = np.int64)
np.logical_and(field[1:], field[:-1], out = above[2:, 1:])
np.cumsum(above, axis = 0, out = above)
np.cumsum(above, axis = 1, out = above)

right = np.zeros((N+1, M+1), dtype = np.int64)
np.logical_and(field[:, 1:], field[:, :-1], out = right[1:, 2:])
np.cumsum(right, axis = 0, out = right)
np.cumsum(right, axis = 1, out = right)

field_csum = np.zeros((N+1, M+1), dtype = np.int64)
field_csum[1:, 1:] = np.cumsum(np.cumsum(field, axis = 1), axis = 0)

for x1, y1, x2, y2 in xyxys:
    x1 -= 1
    y1 -= 1
    ans = field_csum[x2][y2] - field_csum[x1][y2] - field_csum[x2][y1] + field_csum[x1][y1]
    
    y1 += 1
    ans -= right[x2][y2] - right[x1][y2] - right[x2][y1] + right[x1][y1]
    y1 -= 1
    
    x1 += 1
    ans -= above[x2][y2] - above[x1][y2] - above[x2][y1] + above[x1][y1]
    x1 -= 1
    print(ans)

Solution (Python):
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")


n,m,q = list(map(int, input().split()))
s = [None]*n
for i in range(n):
    s[i] = list(map(int, input()))
e1 = [[0]*(m) for _ in range(n-1)]
e2 = [[0]*(m-1) for _ in range(n)]
for i in range(n):
    for j in range(m):
        if i<n-1 and s[i][j]==s[i+1][j]==1:
            e1[i][j] = 1
        if j<m-1 and s[i][j]==s[i][j+1]==1:
            e2[i][j] = 1
def cumsum2(v):
    
    h,w = len(v), len(v[0])
    v.insert(0, [0]*w)
    for l in v:
        l.insert(0, 0)
    for i in range(1,h+1):
        for j in range(1,w+1):
            v[i][j] = v[i-1][j] + v[i][j-1] - v[i-1][j-1] + v[i][j]
    return v
s = cumsum2(s)
if e1:
    e1 = cumsum2(e1)
if e2:
    e2 = cumsum2(e2)
ans = [None]*q
for i in range(q):
    x1,y1,x2,y2 = map(lambda x: int(x)-1, input().split())
    val = s[x2+1][y2+1] - s[x1][y2+1] - s[x2+1][y1] + s[x1][y1]
    if e1:
        val2 = e1[x2][y2+1] - e1[x2][y1] - e1[x1][y2+1] + e1[x1][y1]
    else:
        val2 = 0
    if e2:
        val3 = e2[x2+1][y2] - e2[x2+1][y1] - e2[x1][y2] + e2[x1][y1]
    else:
        val3 = 0
    ans[i] = val - (val2+val3)
write("\n".join(map(str, ans)))

Solution (Python):
import os
import sys

import numpy as np


def solve(n, m, q, sss, xy):
    dp1 = np.zeros((n + 1, m + 1), dtype=np.int64)
    dpv = np.zeros((n + 1, m + 1), dtype=np.int64)
    dph = np.zeros((n + 1, m + 1), dtype=np.int64)
    for i in range(n):
        for j in range(m):
            if sss[i, j] == 1:
                dp1[i + 1, j + 1] = 1

                if i > 0 and sss[i - 1, j] == 1:
                    dpv[i + 1, j + 1] = 1
                if j > 0 and sss[i, j - 1] == 1:
                    dph[i + 1, j + 1] = 1

    for i in range(1, n + 1):
        dp1[i] = dp1[i].cumsum()
        dpv[i] = dpv[i].cumsum()
        dph[i] = dph[i].cumsum()

    for j in range(1, m + 1):
        dp1[:, j] = dp1[:, j].cumsum()
        dpv[:, j] = dpv[:, j].cumsum()
        dph[:, j] = dph[:, j].cumsum()

    buf = np.zeros(q, dtype=np.int64)
    for i in range(q):
        x1 = xy[i * 4 + 0]
        y1 = xy[i * 4 + 1]
        x2 = xy[i * 4 + 2]
        y2 = xy[i * 4 + 3]
        cnt1 = dp1[x2, y2] - dp1[x1 - 1, y2] - dp1[x2, y1 - 1] + dp1[x1 - 1, y1 - 1]
        ver1 = dpv[x2, y2] - dpv[x1, y2] - dpv[x2, y1 - 1] + dpv[x1, y1 - 1]
        hor1 = dph[x2, y2] - dph[x1 - 1, y2] - dph[x2, y1] + dph[x1 - 1, y1]
        buf[i] = cnt1 - ver1 - hor1

    return buf


if sys.argv[-1] == 'ONLINE_JUDGE':
    from numba.pycc import CC

    cc = CC('my_module')
    cc.export('solve', '(i8,i8,i8,i1[:,:],i8[:])')(solve)
    cc.compile()
    exit()

if os.name == 'posix':
    
    from my_module import solve
else:
    from numba import njit

    solve = njit('(i8,i8,i8,i1[:,:],i8[:])', cache=True)(solve)

n, m, q = map(int, input().split())
sss = np.array([list(map(int, input())) for _ in range(n)], dtype=np.int8)
xy = np.fromstring(sys.stdin.read(), dtype=np.int64, sep=' ')
ans = solve(n, m, q, sss, xy)
print('\n'.join(map(str, ans)))

Solution (Python):
import os
import sys

import numpy as np


def solve(n, m, q, sss, xy):
    dp1 = np.zeros((n + 1, m + 1), dtype=np.int64)
    dpv = np.zeros((n + 1, m + 1), dtype=np.int64)
    dph = np.zeros((n + 1, m + 1), dtype=np.int64)
    for i in range(n):
        for j in range(m):
            if sss[i, j] == 1:
                dp1[i + 1, j + 1] = 1

                if i > 0 and sss[i - 1, j] == 1:
                    dpv[i + 1, j + 1] = 1
                if j > 0 and sss[i, j - 1] == 1:
                    dph[i + 1, j + 1] = 1

    for i in range(1, n + 1):
        for j in range(1, m + 1):
            dp1[i, j] += dp1[i, j - 1]
            dpv[i, j] += dpv[i, j - 1]
            dph[i, j] += dph[i, j - 1]

    for j in range(1, m + 1):
        for i in range(1, n + 1):
            dp1[i, j] += dp1[i - 1, j]
            dpv[i, j] += dpv[i - 1, j]
            dph[i, j] += dph[i - 1, j]

    buf = np.zeros(q, dtype=np.int64)
    for i in range(q):
        x1 = xy[i * 4 + 0]
        y1 = xy[i * 4 + 1]
        x2 = xy[i * 4 + 2]
        y2 = xy[i * 4 + 3]
        cnt1 = dp1[x2][y2] - dp1[x1 - 1][y2] - dp1[x2][y1 - 1] + dp1[x1 - 1][y1 - 1]
        ver1 = dpv[x2][y2] - dpv[x1][y2] - dpv[x2][y1 - 1] + dpv[x1][y1 - 1]
        hor1 = dph[x2][y2] - dph[x1 - 1][y2] - dph[x2][y1] + dph[x1 - 1][y1]
        buf[i] = cnt1 - ver1 - hor1

    return buf


if sys.argv[-1] == 'ONLINE_JUDGE':
    from numba.pycc import CC

    cc = CC('my_module')
    cc.export('solve', '(i8,i8,i8,i1[:,:],i8[:])')(solve)
    cc.compile()
    exit()

if os.name == 'posix':
    
    from my_module import solve
else:
    from numba import njit

    solve = njit('(i8,i8,i8,i1[:,:],i8[:])', cache=True)(solve)

n, m, q = map(int, input().split())
sss = np.array([list(map(int, input())) for _ in range(n)], dtype=np.int8)
xy = np.fromstring(sys.stdin.read(), dtype=np.int64, sep=' ')
ans = solve(n, m, q, sss, xy)
print('\n'.join(map(str, ans)))

Solution (Python):
import sys

n, m, q = map(int, input().split())
dp1 = [[0] * (m + 1) for _ in range(n + 1)]
dpv = [[0] * (m + 1) for _ in range(n + 1)]
dph = [[0] * (m + 1) for _ in range(n + 1)]

prev_s = '0' * m

for i in range(n):
    s = input()
    for j in range(m):
        is1 = s[j] == '1'
        additional = 0

        if is1:
            dp1[i + 1][j + 1] = 1

            if i > 0 and prev_s[j] == '1':
                dpv[i + 1][j + 1] = 1
            if j > 0 and s[j - 1] == '1':
                dph[i + 1][j + 1] = 1
    prev_s = s

for i in range(1, n + 1):
    for j in range(1, m + 1):
        dp1[i][j] += dp1[i][j - 1]
        dpv[i][j] += dpv[i][j - 1]
        dph[i][j] += dph[i][j - 1]

for j in range(1, m + 1):
    for i in range(1, n + 1):
        dp1[i][j] += dp1[i - 1][j]
        dpv[i][j] += dpv[i - 1][j]
        dph[i][j] += dph[i - 1][j]

mp = map(int, sys.stdin.read().split())
buf = []
for x1, y1, x2, y2 in zip(mp, mp, mp, mp):
    cnt1 = dp1[x2][y2] - dp1[x1 - 1][y2] - dp1[x2][y1 - 1] + dp1[x1 - 1][y1 - 1]
    ver1 = dpv[x2][y2] - dpv[x1][y2] - dpv[x2][y1 - 1] + dpv[x1][y1 - 1]
    hor1 = dph[x2][y2] - dph[x1 - 1][y2] - dph[x2][y1] + dph[x1 - 1][y1]

    buf.append(cnt1 - ver1 - hor1)

print('\n'.join(map(str, buf)))

Solution (Python):
import sys
input = sys.stdin.readline
N, M, Q = map(int, input().split())
S = [list(map(int, list(input())[: -1])) for _ in range(N)]
e = [[0] * M for _ in range(N + 1)]
for i in range(1, N + 1):
  for j in range(M - 1): e[i][j + 1] = e[i][j] + ((S[i - 1][j] == 1) and (S[i - 1][j + 1] == 1))
for i in range(1, N):
  for j in range(M): e[i + 1][j] += e[i][j]

cs = [[0] * (M + 1) for _ in range(N + 1)]
for i in range(1, N + 1):
  for j in range(M): cs[i][j + 1] = cs[i][j] + (S[i - 1][j] == 1)
for i in range(1, N):
  for j in range(M + 1): cs[i + 1][j] += cs[i][j]

T = [[0] * N for _ in range(M)]
for i in range(N):
  for j in range(M): T[j][i] = S[i][j]

N, M = M, N
ee = [[0] * M for _ in range(N + 1)]
for i in range(1, N + 1):
  for j in range(M - 1): ee[i][j + 1] = ee[i][j] + ((T[i - 1][j] == 1) and (T[i - 1][j + 1] == 1))
for i in range(1, N):
  for j in range(M): ee[i + 1][j] += ee[i][j]




for i in range(Q):
  u, l, d, r = map(lambda x: int(x) - 1, input().split())
  rese = e[d + 1][r] - e[d + 1][l] - e[u][r] + e[u][l]
  resee = ee[r + 1][d] - ee[r + 1][u] - ee[l][d] + ee[l][u]
  res = cs[d + 1][r + 1] - cs[d + 1][l] - cs[u][r + 1] + cs[u][l]
  print(res - rese - resee)

Solution (Python):
mod = 1000000007
eps = 10**-9


def main():
    import sys
    input = sys.stdin.readline

    H, W, Q = map(int, input().split())
    grid = []
    for _ in range(H):
        grid.append(input().rstrip('\n'))

    V = [[0] * (W+1) for _ in range(H+1)]
    for h in range(H):
        for w in range(W):
            V[h+1][w+1] = V[h+1][w] + V[h][w+1] - V[h][w] + int(grid[h][w])
    EH = [[0] * (W+1) for _ in range(H)]
    for h in range(H-1):
        for w in range(W):
            if grid[h][w] == grid[h+1][w] == "1":
                add = 1
            else:
                add = 0
            EH[h+1][w+1] = EH[h+1][w] + EH[h][w+1] - EH[h][w] + add
    EW = [[0] * (W + 1) for _ in range(H + 1)]
    for h in range(H):
        for w in range(W-1):
            if grid[h][w] == grid[h][w+1] == "1":
                add = 1
            else:
                add = 0
            EW[h + 1][w + 1] = EW[h + 1][w] + EW[h][w + 1] - EW[h][w] + add
    for _ in range(Q):
        h1, w1, h2, w2 = map(int, input().split())
        v = V[h2][w2] - V[h2][w1-1] - V[h1-1][w2] + V[h1-1][w1-1]
        eh = EH[h2-1][w2] - EH[h2-1][w1-1] - EH[h1-1][w2] + EH[h1-1][w1-1]
        ew = EW[h2][w2-1] - EW[h2][w1-1] - EW[h1 - 1][w2-1] + EW[h1 - 1][w1-1]
        print(v - eh - ew)


if __name__ == '__main__':
    main()

Solution (Python):
import sys
import numpy as np

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

def precompute(S):
    H, W = S.shape
    V = np.zeros((H + 1, W + 1), np.int64)
    V[1:, 1:] = np.cumsum(S, axis=0).cumsum(axis=1)
    Ex = np.zeros((H, W + 1), np.int64)
    Ex[1:, 1:] = np.cumsum(S[:-1] & S[1:], axis=0).cumsum(axis=1)
    Ey = np.zeros((H + 1, W), np.int64)
    Ey[1:, 1:] = np.cumsum(S[:, :-1] & S[:, 1:], axis=0).cumsum(axis=1)
    return V, Ex, Ey

def main(V, Ex, Ey, query):
    for i in range(len(query) // 4):
        a, b, c, d = query[4 * i:4 * i + 4]
        v = V[c, d] + V[a - 1, b - 1] - V[a - 1, d] - V[c, b - 1]
        e1 = Ex[c - 1, d] + Ex[a - 1, b - 1] - Ex[a - 1, d] - Ex[c - 1, b - 1]
        e2 = Ey[c, d - 1] + Ey[a - 1, b - 1] - Ey[a - 1, d - 1] - Ey[c, b - 1]
        print(v - e1 - e2)

if sys.argv[-1] == 'ONLINE_JUDGE':
    import numba
    from numba.pycc import CC
    i8 = numba.from_dtype(np.int64)
    signature = (i8[:, :], i8[:, :], i8[:, :], i8[:])

    cc = CC('my_module')
    cc.export('main', signature)(main)
    cc.compile()

from my_module import main

H, W, Q = map(int, readline().split())
S = []
for _ in range(H):
    S.append(list(readline().rstrip()))
query = np.array(read().split(), np.int64)

S = (np.concatenate(S) - ord('0')).reshape(H, W)
V, Ex, Ey = precompute(S)

main(V, Ex, Ey, query)

Solution (Python):
n, m, q = map(int, input().split())
s = [input() for _ in range(n)]
a = [[0] * (m + 1) for _ in range(n + 1)]
b = [[0] * (m + 1) for _ in range(n + 1)]
c = [[0] * (m + 1) for _ in range(n + 1)]
for i in range(n):
    for j in range(m):
        a[i + 1][j + 1] = a[i + 1][j]
        if s[i][j] == '1':
            a[i + 1][j + 1] += 1
for i in range(m + 1):
    for j in range(n):
        a[j + 1][i] += a[j][i]
for i in range(n):
    for j in range(1, m):
        b[i + 1][j + 1] = b[i + 1][j]
        if s[i][j] == s[i][j - 1] == '1':
            b[i + 1][j + 1] += 1
for i in range(m + 1):
    for j in range(n):
        b[j + 1][i] += b[j][i]
for j in range(m):
    for i in range(1, n):
        c[i + 1][j + 1] = c[i][j + 1]
        if s[i][j] == s[i - 1][j] == '1':
            c[i + 1][j + 1] += 1
for j in range(n + 1):
    for i in range(m):
        c[j][i + 1] += c[j][i]
for _ in range(q):
    x1, y1, x2, y2 = map(int, input().split())
    x1 -= 1
    y1 -= 1
    print((a[x2][y2] - a[x1][y2] - a[x2][y1] + a[x1][y1]) -
          (b[x2][y2] - b[x1][y2] - b[x2][y1 + 1] + b[x1][y1 + 1]) -
          (c[x2][y2] - c[x1 + 1][y2] - c[x2][y1] + c[x1 + 1][y1]))

Solution (Python):
def examA():
    N, A, B = LI()
    if N==1:
        if A==B:
            print(1)
        else:
            print(0)
        return
    ans = (N-2)*max(B-A,0) + 1
    if B<A:
        ans = 0
    print(ans)
    return

def examB():
    S = SI()
    N = len(S)
    cnt = 0
    for i,s in enumerate(S):
        cnt += N-1
        if s=="U":
            cnt += i
        else:
            cnt += (N-i-1)
    ans = cnt
    print(ans)
    return

def examC():
    N, M, q = LI()
    S = [SI()for _ in range(N)]
    Q = [LI()for _ in range(q)]
    S1 = [[0]*(M+1) for _ in range(N+1)]
    S2 = [[0]*(M+1) for _ in range(N+1)]
    S3 = [[0]*(M+1) for _ in range(N+1)]
    for i in range(N):
        now = 0
        for j in range(M):
            now += int(S[i][j])
            S1[i+1][j+1] = S1[i][j+1]+now
    for i in range(N):
        now = 0
        for j in range(M-1):
            if S[i][j]=="1" and S[i][j+1]=="1":
                now += 1
            S2[i+1][j+1] += S2[i][j+1]+now
    for i in range(M):
        now = 0
        for j in range(N-1):
            if S[j][i]=="1" and S[j+1][i]=="1":
                now += 1
            S3[j+1][i+1] += S3[j+1][i]+now
    ans = [0]*q
    
    
    for i in range(q):
        x1, y1, x2, y2 = Q[i]
        
        x1 -= 1; y1 -= 1
        cur = S1[x2][y2] - S1[x1][y2] - S1[x2][y1] + S1[x1][y1]
        cur -= S2[x2][y2-1] - S2[x1][y2-1] - S2[x2][y1] + S2[x1][y1]
        cur -= S3[x2-1][y2] - S3[x1][y2] - S3[x2-1][y1] + S3[x1][y1]
        ans[i] = cur
    for v in ans:
        print(v)
    return

def examD():
    A = I()
    B = I()

    ans = 0
    print(ans)
    return

def examE():
    ans = 0
    print(ans)
    return

def examF():
    ans = 0
    print(ans)
    return

def test():

    return

import sys,bisect,itertools,heapq,math,random
from copy import deepcopy
from heapq import heappop,heappush,heapify
from collections import Counter,defaultdict,deque
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
def I(): return int(readline())
def LI(): return list(map(int,sys.stdin.readline().split()))
def LSI(): return list(map(str,sys.stdin.readline().split()))
def LS(): return sys.stdin.readline().split()
def SI(): return sys.stdin.readline().strip()
global mod,mod2,inf,alphabet,_ep
mod = 10**9 + 7
mod2 = 998244353
inf = 10**18
_ep = 10**(-12)
alphabet = [chr(ord('a') + i) for i in range(26)]

sys.setrecursionlimit(10**7)

if __name__ == '__main__':
    examC()

"""
3 4 4 
1101 
0110 
1101 
1 1 3 4 
1 1 3 1 
2 2 3 4 
1 2 2 4 
"""

Solution (Python):
import sys
input = sys.stdin.readline

n, m, q = map(int, input().split())
S = [[0]*(m+1)] + [[0] + list(map(int, list(input().rstrip()))) for _ in range(n)]
B = [[0]*(m+1) for _ in range(n+1)]
U = [[0]*(m+1) for _ in range(n+1)]
L = [[0]*(m+1) for _ in range(n+1)]
for i in range(1, n+1):
  for j in range(1, m+1):
    B[i][j] = B[i-1][j] + B[i][j-1] - B[i-1][j-1] + S[i][j]
    U[i][j] = U[i-1][j] + U[i][j-1] - U[i-1][j-1] + U[i][j]
    L[i][j] = L[i-1][j] + L[i][j-1] - L[i-1][j-1] + L[i][j]
    if S[i][j]:
      U[i][j] += S[i-1][j]
      L[i][j] += S[i][j-1]
    
def calc_b(x1, y1, x2, y2):
  return B[x2][y2] - B[x1-1][y2] - B[x2][y1-1] + B[x1-1][y1-1]
def calc_u(x1, y1, x2, y2):
  return U[x2][y2] - U[x1][y2] - U[x2][y1-1] + U[x1][y1-1]
def calc_l(x1, y1, x2, y2):
  return L[x2][y2] - L[x1-1][y2] - L[x2][y1] + L[x1-1][y1]

for _ in range(q):
  x1, y1, x2, y2 = map(int, input().split())
  res = calc_b(x1, y1, x2, y2) - (calc_u(x1, y1, x2, y2) + calc_l(x1, y1, x2, y2))
  print(res)

Solution (Python):
import sys

input=sys.stdin.readline

N,M,Q=map(int,input().split())
S=[]
for i in range(N):
    S.append(input())

blue=[[0 for j in range(M+1)] for i in range(N+1)]
for i in range(1,N+1):
    for j in range(1,M+1):
        blue[i][j]+=blue[i][j-1]+int(S[i-1][j-1])

for j in range(1,M+1):
    for i in range(1,N+1):
        blue[i][j]+=blue[i-1][j]

yoko=[[0 for j in range(M+1)] for i in range(N+1)]
for i in range(1,N+1):
    for j in range(2,M+1):
        yoko[i][j]+=int(S[i-1][j-1])*int(S[i-1][j-2])+yoko[i][j-1]

for j in range(M+1):
    for i in range(1,N+1):
        yoko[i][j]+=yoko[i-1][j]

tate=[[0 for i in range(M+1)] for j in range(N+1)]
for j in range(1,M+1):
    for i in range(2,N+1):
        tate[i][j]+=tate[i-1][j]+int(S[i-2][j-1])*int(S[i-1][j-1])

for i in range(N+1):
    for j in range(1,M+1):
        tate[i][j]+=tate[i][j-1]

for _ in range(Q):
    x1,y1,x2,y2=map(int,input().split())
    res=(blue[x2][y2]-blue[x1-1][y2])-(blue[x2][y1-1]-blue[x1-1][y1-1])
    res-=(tate[x2][y2]-tate[x2][y1-1])-(tate[x1][y2]-tate[x1][y1-1])
    res-=(yoko[x2][y2]-yoko[x1-1][y2])-(yoko[x2][y1]-yoko[x1-1][y1])
    print(res)

Solution (Python):
def main():
    from itertools import accumulate as ac
    n, m, q = map(int, input().split())
    ten = [list(map(lambda x:x == "1", list(input()))) for _ in [0]*n]
    xy = [list(map(int, input().split())) for _ in [0]*q]
    ans = [0]*q

    hen = [[0]*m]+[None]*n
    for i in range(n):
        temp = [0]
        for j in range(m-1):
            temp.append(ten[i][j] & ten[i][j+1])
        hen[i+1] = list(ac(temp))
    for i in range(n):
        for j in range(1, m):
            hen[i+1][j] += hen[i][j]

    for i, (x1, y1, x2, y2) in enumerate(xy):
        ans[i] = -hen[x2][y2-1]-hen[x1-1][y1-1] + \
            hen[x2][y1-1]+hen[x1-1][y2-1]

    hen = [[0]*(m+1)]+[None]*(n-1)
    for i in range(n-1):
        temp = [0]
        for j in range(m):
            temp.append(ten[i][j] & ten[i+1][j])
        hen[i+1] = list(ac(temp))
    for i in range(n-1):
        for j in range(1, m+1):
            hen[i+1][j] += hen[i][j]

    for i, (x1, y1, x2, y2) in enumerate(xy):
        ans[i] -= hen[x2-1][y2]+hen[x1-1][y1-1] - \
            hen[x2-1][y1-1]-hen[x1-1][y2]

    del hen

    ten = [[0]*(m+1)]+ten
    for i in range(1, n+1):
        ten[i] = [0]+list(ac(ten[i]))
    for i in range(n):
        for j in range(1, m+1):
            ten[i+1][j] += ten[i][j]

    for i, (x1, y1, x2, y2) in enumerate(xy):
        print(ans[i]+ten[x2][y2]+ten[x1-1][y1-1]-ten[x2][y1-1]-ten[x1-1][y2])


main()

Solution (Python):
import gc


class Ruiseki2D():
    """二次元累積和を構築する
    構築O(HW), 区間取得クエリO(1)
    """
    def __init__(self, matrix):
        self.h = len(matrix)
        self.w = len(matrix[0])
        self.ruiseki = [[0] * (self.w + 1) for _ in range(self.h + 1)]

        for i in range(self.h):
            for j in range(self.w):
                self.ruiseki[i + 1][j + 1] = self.ruiseki[i + 1][j] + matrix[i][j]
        for i in range(self.h):
            for j in range(self.w):
                self.ruiseki[i + 1][j + 1] += self.ruiseki[i][j + 1]

    def get_sum(self, hl, hr, wl, wr):
        """[hl, hr), [wl, wr) で囲まれた部分の和を求める"""
        return self.ruiseki[hr][wr] + self.ruiseki[hl][wl] - self.ruiseki[hr][wl] - self.ruiseki[hl][wr]


n, m, q = map(int, input().split())
s = [list(map(int, input())) for i in range(n)]
query = [list(map(int, input().split())) for i in range(q)]

ru_col_cnt = Ruiseki2D(s)


yoko_edge = [[0] * (m + 1) for i in range(n + 1)]
for i in range(n):
    for j in range(m - 1):
        if s[i][j] == 1 and s[i][j + 1] == 1:
            yoko_edge[i + 1][j + 1] = 1
for i in range(n):
    for j in range(m):
        yoko_edge[i + 1][j + 1] += yoko_edge[i + 1][j]
for i in range(n):
    for j in range(m):
        yoko_edge[i + 1][j + 1] += yoko_edge[i][j + 1]

tate_edge = [[0] * (m + 1) for i in range(n + 1)]
for i in range(n - 1):
    for j in range(m):
        if s[i][j] == 1 and s[i + 1][j] == 1:
            tate_edge[i + 1][j + 1] = 1
for i in range(n):
    for j in range(m):
        tate_edge[i + 1][j + 1] += tate_edge[i + 1][j]
for i in range(n):
    for j in range(m):
        tate_edge[i + 1][j + 1] += tate_edge[i][j + 1]



def yoko_get_sum(hl, hr, wl, wr):
    return yoko_edge[hr][wr] + yoko_edge[hl][wl] - yoko_edge[hr][wl] - yoko_edge[hl][wr]
  
def tate_get_sum(hl, hr, wl, wr):
    return tate_edge[hr][wr] + tate_edge[hl][wl] - tate_edge[hr][wl] - tate_edge[hl][wr]

for hl, wl, hr, wr in query:
    hl, wl = hl - 1, wl - 1
    ans = ru_col_cnt.get_sum(hl, hr, wl, wr)
    ans -= yoko_get_sum(hl, hr, wl, wr - 1) + tate_get_sum(hl, hr - 1, wl, wr)
    print(ans)

Solution (Python):
def get(a, b, c, d, dp):
    return dp[c][d]-dp[a][d]-dp[c][b]+dp[a][b]

n, m, q = map(int, input().split())
s = [list(map(int, list(input()))) for _ in range(n)]
vers = [[0]*(m+1) for i in range(n+1)]
edges = [[0]*(m+1) for i in range(n+1)]
tate = [[0]*(m+1) for i in range(n+1)]
yoko = [[0]*(m+1) for i in range(n+1)]

for i in range(n):
    for j in range(m):
        if s[i][j] == 0:
            continue
        vers[i+1][j+1] = 1
        nx, ny = i+1, j
        if 0<=nx<n and 0<=ny<m and s[nx][ny] == 1:
            edges[i+1][j+1]+=1
            edges[i+2][j+1]+=1
            tate[i+1][j+1]+=1
        nx, ny = i, j+1
        if 0<=nx<n and 0<=ny<m and s[nx][ny] == 1:
            edges[i+1][j+1]+=1
            edges[i+1][j+2]+=1
            yoko[i+1][j+1]+=1


for i in range(1, n+1):
    for j in range(1, m+1):
            vers[i][j]+=vers[i][j-1]
            edges[i][j]+=edges[i][j-1]
            tate[i][j]+=tate[i][j-1]
for j in range(1, m+1):
    for i in range(1, n+1):
            vers[i][j]+=vers[i-1][j]
            edges[i][j]+=edges[i-1][j]
            yoko[i][j]+=yoko[i-1][j]
for _ in range(q):
    a, b, c, d = map(int, input().split())
    a-=1
    b-=1
    vcnt = get(a, b, c, d, vers)
    ecnt = get(a, b, c, d, edges)
    ecnt-=(yoko[c][b]-yoko[a][b])
    ecnt-=(tate[a][d]-tate[a][b])
    ecnt-=(yoko[c][d]-yoko[a][d])
    ecnt-=(tate[c][d]-tate[c][b])
    ecnt//=2
    print(vcnt-ecnt)

Solution (Python):
def solve():
    import sys
    input = sys.stdin.readline

    N, M, Q = map(int, input().split())
    Sss = [input().rstrip() for _ in range(N)]

    numVertexss = [[0]*(M+1)] + [[0] + list(map(int, Ss)) for Ss in Sss]
    numEdgeVss = [[0]*(M+1) for _ in range(N+1)]
    numEdgeHss = [[0]*(M+1) for _ in range(N+1)]
    for i in range(1, N):
        for j in range(1, M+1):
            if numVertexss[i][j] and numVertexss[i+1][j]:
                numEdgeVss[i][j] = 1
    for i in range(1, N+1):
        for j in range(1, M):
            if numVertexss[i][j] and numVertexss[i][j+1]:
                numEdgeHss[i][j] = 1

    def getAccAss(Ass):
        for x in range(1, N+1):
            for y in range(M+1):
                Ass[x][y] += Ass[x-1][y]
        for x in range(N+1):
            for y in range(1, M+1):
                Ass[x][y] += Ass[x][y-1]
    def getRangeSum2D(accAss, xFr, xTo, yFr, yTo):
        return accAss[xTo+1][yTo+1] - accAss[xTo+1][yFr] - accAss[xFr][yTo+1] + accAss[xFr][yFr]

    getAccAss(numVertexss)
    getAccAss(numEdgeVss)
    getAccAss(numEdgeHss)

    anss = []
    for _ in range(Q):
        x1, y1, x2, y2 = map(int, input().split())
        x1, y1, x2, y2 = x1-1, y1-1, x2-1, y2-1
        ans = getRangeSum2D(numVertexss, x1, x2, y1, y2)
        ans -= getRangeSum2D(numEdgeVss, x1, x2-1, y1, y2)
        ans -= getRangeSum2D(numEdgeHss, x1, x2, y1, y2-1)
        anss.append(ans)

    print('\n'.join(map(str, anss)))


solve()

Solution (Python):
def main():
    import itertools
    import os
    import sys

    if os.getenv("LOCAL"):
        sys.stdin = open("_in.txt", "r")

    sys.setrecursionlimit(10 ** 9)
    INF = float("inf")
    IINF = 10 ** 18
    MOD = 10 ** 9 + 7


    


    def cumsum(it):
        """
        累積和
        :param collections.Iterable it:
        """
        cs = 0
        ret = [0]
        for v in it:
            cs += v
            ret.append(cs)
        return ret


    def cumsum_mat(mat):
        h = len(mat)
        w = len(mat[0])

        tmp = []
        for row in mat:
            tmp.append(cumsum(row))
        del mat

        ret = [[0] * (w + 1) for _ in range(h + 1)]
        for c in range(w + 1):
            col = [row[c] for row in tmp]
            for r, a in enumerate(cumsum(col)):
                ret[r][c] = a
        return ret


    N, M, Q = list(map(int, sys.stdin.buffer.readline().split()))
    S = [list(map(int, sys.stdin.buffer.readline().decode().rstrip())) for _ in range(N)]
    XY = [list(map(int, sys.stdin.buffer.readline().split())) for _ in range(Q)]


    def count_v(x1, y1, x2, y2):
        return counts_cum[x2][y2] - counts_cum[x1 - 1][y2] - counts_cum[x2][y1 - 1] + counts_cum[x1 - 1][y1 - 1]


    def count_ev(x1, y1, x2, y2):
        return edges_v_cum[x2][y2] - edges_v_cum[x1][y2] - edges_v_cum[x2][y1 - 1] + edges_v_cum[x1][y1 - 1]


    def count_eh(x1, y1, x2, y2):
        return edges_h_cum[x2][y2] - edges_h_cum[x1 - 1][y2] - edges_h_cum[x2][y1] + edges_h_cum[x1 - 1][y1]


    vs = []
    e1 = []
    e2 = []

    counts_cum = cumsum_mat(S)
    for x1, y1, x2, y2 in XY:
        vs.append(count_v(x1, y1, x2, y2))
    del counts_cum

    edges_v = [[0] * M for _ in range(N)]
    for h, w in itertools.product(range(1, N), range(M)):
        edges_v[h][w] = S[h][w] & S[h - 1][w]
    edges_v_cum = cumsum_mat(edges_v)
    for x1, y1, x2, y2 in XY:
        e1.append(count_ev(x1, y1, x2, y2))
    del edges_v
    del edges_v_cum

    edges_h = [[0] * M for _ in range(N)]
    for h, w in itertools.product(range(N), range(1, M)):
        edges_h[h][w] = S[h][w] & S[h][w - 1]
    edges_h_cum = cumsum_mat(edges_h)
    for x1, y1, x2, y2 in XY:
        e2.append(count_eh(x1, y1, x2, y2))
    del edges_h
    del edges_h_cum

    ans = [v - e1 - e2 for v, e1, e2 in zip(vs, e1, e2)]
    print(*ans, sep='\n')


main()

Solution (Python):
import sys
input = sys.stdin.readline

N, M, Q = map(int, input().split())
state = [[int(a) for a in list(input().rstrip())] for _ in range(N)]
Query = [list(map(int, input().split())) for _ in range(Q)]
updown = [[0]*(M+1) for _ in range(N)]
for h in range(N-1):
    st = 0
    for w in range(M):
        if state[h][w] and state[h+1][w]:
            st += 1
        updown[h+1][w+1] = st

rileft = [[0]*M for _ in range(N+1)]
for w in range(M-1):
    st = 0
    for h in range(N):
        if state[h][w] and state[h][w+1]:
            st += 1
        rileft[h+1][w+1] = st

dp = [[0]*(M+1) for _ in range(N+1)]
for h in range(N):
    have = False
    connected = 0
    for w in range(M):
        if state[h][w]:
            if not have:
                connected += 1
                have = True
        else:
            if have:
                have = False
        dp[h+1][w+1] = dp[h][w+1] + connected - updown[h][w+1]

for n1, m1, n2, m2 in Query:
    midcon = dp[n2][m2] - dp[n1-1][m2] + updown[n1-1][m2]
    lecon = dp[n2][m1-1] - dp[n1-1][m1-1] + updown[n1-1][m1-1]
    ans = midcon - lecon + (rileft[n2][m1-1] - rileft[n1-1][m1-1])
    print(ans)

Solution (Python):
import sys
from itertools import accumulate
readline = sys.stdin.readline

def accumulate2d(X):
    N = len(X)
    M = len(X[0])
    
    for i in range(1, N):
        for j in range(1, M):
            X[i][j] += X[i][j-1]
    
    for j in range(1, M):
        for i in range(1, N):
            X[i][j] += X[i-1][j]
    
    return X
    

N, M, Q = map(int, readline().split())
G = [[0]*(M+2)] + [[0] + list(map(int, readline().strip())) + [0] for _ in range(N)] + [[0]*(M+2)]
U = [[0]*(M+2) for _ in range(N+2)]
L = [[0]*(M+2) for _ in range(N+2)]
for i in range(1, N+1):
    for j in range(1, M+1):
        if G[i][j]:
            if G[i-1][j]:
                U[i][j] += 1
            if G[i][j-1]:
                L[i][j] += 1
G = accumulate2d(G)
L = accumulate2d(L)
U = accumulate2d(U)
Ans = [None]*Q 
for qu in range(Q):
    x1, y1, x2, y2 = map(int, readline().split())
    u = U[x1][y1-1] + U[x2][y2] - U[x1][y2] - U[x2][y1-1]
    l = L[x1-1][y1] + L[x2][y2] - L[x1-1][y2] - L[x2][y1]
    g = G[x1-1][y1-1] + G[x2][y2] - G[x1-1][y2] - G[x2][y1-1]
    Ans[qu] = g - u - l
print('\n'.join(map(str, Ans)))