x=(r[j]*(r[j]+1)/2-l[j]*(l[j]-1)/2)*(r[i]-l[i]+1);//(sum of numbers from l[j] to r[j])*(number of calumn in a given submatrix)

y=(r[i]*(r[i]+1)/2-l[i]*(l[i]-1)/2)*(r[j]-l[j]+1);

cout<<x+y%k;

The following idea may help in solving the problem:

Suppose that the function T( n, m, p ) represents the sum of all elements in a sub-matrix with n rows and m columns whose element in row i and column j is equal to (p+i+j-2) mod K, and the summation is computed modulo R = 1000000007. The first element of this sub-matrix at i = j = 1 is p mod K, and all elements in the sub-matrix are in the integer number interval U = [0,K-1].

It follows that:

a) For any value of p in U,

T( n, m, p ) = 0, if m < 1 or n < 1

T( K, K, p ) = K K (K-1) / 2, as each row and column of the modulo sub-matrix in this case is a shift-and-rotate version of the vector [ 0 1 2 .... K-1 ] by a distinct shift value u in U.

b) The answer to the query Q( Li, Lj, Ri, Rj ) is equal to T( Ri-Li+1, Rj-Lj+1, (Li+Lj) mod K )

c) If n or m is larger than K-1, then rule (a) can be used to compute T( n, m, p ) in terms of its constituent sub-matrices as follows.

T( n, m, p ) = ( m v+n w-K v w ) K ( K-1 ) / 2 + T( n mod K, m mod K, p ), where v = n div K and w = m div K.

d) T( n, m, 0 ) when n and m are in U can be computed in constant-time using a close-form expression.

The first row and the first column of the sub-matrix are [ 0 1 2 .... m-1 ] and [ 0 1 2 .... n-1 ], respectively. Deriving the closed-form expression is a mathematical exercise.

e) T( n, m, p ) can be computed recursively when p > 0, where n and m are both in U, as the sum of four sub-matrices A, B, C, and D as shown when both n and m are larger than r, where r = K-p:

A B

C D

The sub-matrix A is symmetric with r rows, r columns, and its first row is [ p p+1 p+2 .... K-1 ]. Computing T( r, r, p ) in constant-time can be done using closed-form expression. Again, deriving the closed-form expression is a mathematical exercise. It can be shown that

T( r, r, p ) = r [ K ( r+1 ) / 2-r ]

The sub-matrices B and C have their first elements equal to 0, and the total summation of both sub-matrices is equal to T( r, m-r, 0 ) + T( n-r, r, 0 ). The sub-matrix D has n-r rows m-r columns, its first element is r, and its summation is T( n-r, m-r, r ).

If both n and m are less than or equal to r, then only the sub-matrix A exists, and it is asymmetric when n is not equal to m. If either n or m is less than or equal to r, and the other dimension is larger than r, then only two sub-matrices exist, either (A and B) or (A and C).

f) When m and n are in U, T( m, n, p ) can be computed only once to find the answer to a query, and its value is stored to be used in finding the answer to subsequent queries if needed. All the summation operations are computed modulo R.

Finally, it should be noted that the tuple < m, n, p > is a signature for the modulo sub-matrix as the sum of its elements T( m, n, p ) is invariant to its location in the overall matrix.

Submatrix Sum Queries

Well, I have invented this task :P

You have solution here (on Serbian) : https://algora.petlja.org/t/resenja-zadataka-okruzno-2018/1354/9

My idea was changing big submatrix to submatrix with sides (R_i - L_i + 1)mod K and (R_j - L_j + 1)mod K. It can be easily done because each line (verical or horizontal) with K elements has constant sum (all elements 0, 1..., K - 1 shifted for some amount of places). After it, "smaller" submatrix is easier for summing ( you can divide this smaller submatrix on three parts : each diagonal is larger for one than previous, all diagonals are equal length — has length same as smaller side of "smaller" submatrix, and third part diagonals decreasing ). This is easier for summing because all elements at fixed diagonal are equal and at most once will happen "reset" (diagonals will start again from 0) , you must be careful with it.

Here is code

There is one easier solution for implementation and smaller amount of math. Something similar to sbakic link, but with O(1) memory and one function for calculating value in such sumatrix with upper corner (1, 1). Similar idea for getting "smaller" submatrix, than you have special "smaller" submatrix, with value in upper corner equal to 2. You can sum values over rows, first it will be arithmetic progression of rows with positive coefficient and after it every next row will decrease for constant value than previous.

I hope this was helpful.

The following is a C++17 implementation of the idea suggested yesterday.

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

inline ll& mod( ll &x ) { return x %= 1000000007LL; }

int main()
{
    ll K, L, N, M, Q;

    ios_base::sync_with_stdio( false ), cin.tie( nullptr ), cout.tie( nullptr );

    cin >> N >> M >> K >> Q, L = K, L *= K - 1, L /= 2, mod( L );

    struct modulo_matrix_t
    {
        ll K; modulo_matrix_t( ll k ) : K( k ) {}

        ll T( ll n, ll m, ll p )
        {
            if ( n > m )
            {
                ll q = n; n = m, m = q;
            }

            if ( n < 1 )
                return 0;

            if ( n == 1 && m == 1 )
                return p;

            ll r, u, v, t = K;

            if ( p == 0 )
                t = n - 1, t /= 2,
                mod( t ) *= 2 * n + 1,
                mod( t ) += T( n, m - n, n );
            else
                r = K - p, u = m - r, v = n - r, t *= r + 1, t -= 2 * r,
                mod( t ) *= r, t /= 2,
                mod( t ) += T( r, u, 0 ),
                mod( t ) += T( v, r, 0 ),
                mod( t ) += T( v, u, r );

            return mod( t );
        }

    } modulo_matrix( K );

    for( ll q = 0; q < Q; q++ )
    {
        ll Li, Lj, Ri, Rj; cin >> Li >> Lj >> Ri >> Rj;

        ll n = Ri - Li + 1, m = Rj - Lj + 1, p = ( Li + Lj ) % K;

        ll v = n / K, w = m / K, z = v * w, t = m * v - K * mod( z );

        mod( t ) += n * w, 
        mod( t ) *= L, 
        mod( t ) += modulo_matrix.T( m % K, n % K, p );

        cout << mod( t ) << '\n';
    }
}