We will find the values of \(f(i,j)\) for all \(N^2\) cells. Let \(C_{i,j}\) be the number of participants living in cell \((i,j)\).

Contribution from cell \((i,j)\) to \(f(i+x,j+y)\)

The contribution of \(C_{i,j}\) to \(f(i+x,\ j+y)\), denoted as \(d[x][y]\), can be written as follows:

• \(d[x][y] = \min(\max(|x|,\ |y|),\ N) \times C_{i,j}\)

Now, consider a two-dimensional array \(d'\) whose cumulative sums yield \(d\). In other words, consider the \(d'\) such that \(\displaystyle d[x][y] = NC_{i,j} + \sum_{x' \leq x,\ y' \leq y} d'[x'][y']\) for all \((x,y)\). This can be derived from \(d'[x][y] = d[x][y] - d[x-1][y] - d[x][y-1] + d[x-1][y-1]\), specifically:

• \(d'[k][k] = (-1) \times C_{i,j}\) for all integers \(k\) with \(-N+1 \leq k \leq N\).
• \(d'[k][-k+1] = (+1) \times C_{i,j}\) for all integers \(k\) with \(-N+1 \leq k \leq N\).
• \(d'[x][y] = 0\) for all other places.

This can be expressed as additions onto a constant number of elements, via a special cumulative-sum trick applied in two diagonal directions.

Contribution to all cells.

The contributions of all \(N^2\) cells can be computed by adding the contribution \(d[x][y]\) of cell \((i,j)\) onto the position \((i+x, j+y)\). Therefore, all values of \(f(i,j)\) can be computed by taking the cumulative sums of a two-dimensional array obtained by appropriately distributing the contribution of \(d'\) from each cell, and then globally adding \(N \sum_{i,j} C_{i,j}\).

Sample code (C++)

#include <iostream>
using std::cin;
using std::cout;
using std::cerr;
using std::endl;
#include <vector>
using std::vector;
using std::pair;
#include <map>
using std::map;
using std::max;
using std::min;


#ifdef DEBUG
const int debug = 1;
#else
const int debug = 0;
#endif

using ll = int64_t;
using P = pair<ll, ll>;

const ll FOD = 998244353;
#include <atcoder/modint>
using mint = atcoder::modint998244353;
#include <atcoder/convolution>
using atcoder::convolution;

ll n, k, h, w, m;
vector<ll> a, b;

inline void output (const ll x) {
	cout << x << "\n";
}

void solve() {
	ll c[n][n];
	for (ll i = 0; i < n; i++) {
		for (ll j = 0; j < n; j++) {
			c[i][j] = a[i] * b[j] % m;
		}
	}

	const ll nn = 3*n+1;
	const ll offset = n;
	ll sx[nn][nn], sy[nn][nn], s[nn][nn];
	for (ll i = 0; i < nn; i++) {
		for (ll j = 0; j < nn; j++) {
			sx[i][j] = 0;
			sy[i][j] = 0;
			s[i][j] = 0;
		}
	}

	ll outsum = 0;
	for (ll i = 0; i < n; i++) {
		for (ll j = 0; j < n; j++) {
			outsum += c[i][j] * n;

			sx[(offset+i) - (n-1) - 0][(offset+j) - (n-1) - 0] += c[i][j] * (-1);
			sx[(offset+i) + (n+0) + 1][(offset+j) + (n+0) + 1] -= c[i][j] * (-1);

			sy[(offset+i) - (n-1) - 0][(offset+j) + (n+0) + 0] += c[i][j] * (+1);
			sy[(offset+i) + (n+0) + 1][(offset+j) - (n-1) - 1] -= c[i][j] * (+1);

		}
	}

	// sx, sy
	for (ll i = 0; i < nn; i++) {
		for (ll j = 0; j < nn; j++) {
			if (i - 1 >= 0 && j - 1 >= 0) sx[i][j] += sx[i - 1][j - 1];
			if (i - 1 >= 0 && j + 1 < nn) sy[i][j] += sy[i - 1][j + 1];

			s[i][j] = sx[i][j] + sy[i][j];
		}
	}

	for (ll i = 0; i < nn; i++) {
		ll x = 0;
		for (ll j = 0; j < nn; j++) {
			s[i][j] += x;
			x = s[i][j];
		}
	}
	for (ll j = 0; j < nn; j++) {
		ll x = 0;
		for (ll i = 0; i < nn; i++) {
			s[i][j] += x;
			x = s[i][j];
		}
	}

	ll ans = 0;
	for (ll i = 0; i < n; i++) {
		for (ll j = 0; j < n; j++) {
			ll item = outsum + s[offset+i][offset+j];
			if (debug) {
				cerr << i << " " << j << ": " << item << endl;
			}

			ans ^= (item + i*n + j);
		}
	}

	output(ans);


	return;
}

int main (void) {
	std::cin.tie(nullptr);
	std::ios_base::sync_with_stdio(false);

	cin >> n;
	cin >> m;
	a.resize(n);
	b.resize(n);
	for (ll i = 0; i < n; i++) {
		cin >> a[i];
	}
	for (ll i = 0; i < n; i++) {
		cin >> b[i];
	}

	
	solve();

	return 0;
}