Problem 453: Lattice Quadrilaterals

View on Project Euler

Project Euler Problem 453 Solution

EulerSolve provides an optimized solution for Project Euler Problem 453, Lattice Quadrilaterals, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $$G_{m,n}=\{0,1,\dots,m\}\times\{0,1,\dots,n\},\qquad N=(m+1)(n+1).$$ We seek \(Q(m,n)\), the number of simple quadrilaterals whose four vertices are distinct lattice points of \(G_{m,n}\). The checkpoint values used by the implementations are $$Q(2,2)=94,\quad Q(3,7)=39590,\quad Q(12,3)=309000,\quad Q(123,45)=70542215894646.$$ The actual target is \(Q(12345,6789)\bmod 135707531\). Mathematical Approach The solution does not enumerate quadruples directly. Instead it decomposes the problem into four global quantities: collinear triples, collinear quadruples, the total boundary-point count over all lattice triangles, and the total doubled-area count over all lattice triangles. Step 1: Start from all 4-point sets There are \(\binom{N}{4}\) ways to choose four lattice points. A set is immediately invalid if at least three of its points are collinear, because then no simple quadrilateral can use all four chosen vertices. For a segment with displacement \((i,j)\), the number of lattice points on the segment equals \(\gcd(i,j)+1\). Therefore the number of interior lattice points on that segment is \(\gcd(i,j)-1\). This gcd rule is the basic ingredient behind the whole counting argument. Step 2: Count collinear triples and quadruples by displacement class Let \(\mathcal{D}\) be the set of all nonzero displacements \((i,j)\) with \(0\le i\le m\) and \(0\le j\le n\)....

Detailed mathematical approach

Problem Summary

Let

$$G_{m,n}=\{0,1,\dots,m\}\times\{0,1,\dots,n\},\qquad N=(m+1)(n+1).$$

We seek \(Q(m,n)\), the number of simple quadrilaterals whose four vertices are distinct lattice points of \(G_{m,n}\). The checkpoint values used by the implementations are

$$Q(2,2)=94,\quad Q(3,7)=39590,\quad Q(12,3)=309000,\quad Q(123,45)=70542215894646.$$

The actual target is \(Q(12345,6789)\bmod 135707531\).

Mathematical Approach

The solution does not enumerate quadruples directly. Instead it decomposes the problem into four global quantities: collinear triples, collinear quadruples, the total boundary-point count over all lattice triangles, and the total doubled-area count over all lattice triangles.

Step 1: Start from all 4-point sets

There are \(\binom{N}{4}\) ways to choose four lattice points. A set is immediately invalid if at least three of its points are collinear, because then no simple quadrilateral can use all four chosen vertices.

For a segment with displacement \((i,j)\), the number of lattice points on the segment equals \(\gcd(i,j)+1\). Therefore the number of interior lattice points on that segment is \(\gcd(i,j)-1\). This gcd rule is the basic ingredient behind the whole counting argument.

Step 2: Count collinear triples and quadruples by displacement class

Let \(\mathcal{D}\) be the set of all nonzero displacements \((i,j)\) with \(0\le i\le m\) and \(0\le j\le n\). For \((i,j)\in\mathcal{D}\), define

$$g_{i,j}=\gcd(i,j),$$

with the conventions \(\gcd(i,0)=i\) and \(\gcd(0,j)=j\). The number of segments having that displacement is

$$M_{i,j}=\begin{cases} (n+1)(m+1-i), & j=0,\\ (m+1)(n+1-j), & i=0,\\ 2(m+1-i)(n+1-j), & i,j\ge 1. \end{cases}$$

The factor \(2\) in the last case accounts for the two slopes \((i,j)\) and \((i,-j)\). Then

$$L_3=\sum_{(i,j)\in\mathcal{D}} M_{i,j}\bigl(g_{i,j}-1\bigr),\qquad L_4=\sum_{(i,j)\in\mathcal{D}} M_{i,j}\binom{g_{i,j}-1}{2}.$$

Each interior lattice point on a segment produces one collinear triple with the two endpoints, so \(L_3\) is exactly the number of collinear triples. Choosing two interior points on the same segment produces one collinear quadruple, so \(L_4\) is exactly the number of collinear quadruples.

If every collinear triple is extended by an arbitrary fourth point, we count \((N-3)L_3\) bad 4-point sets. A set of four collinear points contains four different triples, so it has been counted four times and must be corrected by subtracting \(3L_4\). Hence

$$C_{\mathrm{bad}}=(N-3)L_3-3L_4.$$

Step 3: Why triangle interior points matter

Now consider only 4-point sets with no three collinear. Such a set has exactly two possibilities:

1. The four points are in convex position, and there is exactly one simple quadrilateral.

2. One point lies strictly inside the triangle formed by the other three, and there are exactly three simple quadrilaterals.

So every nondegenerate 4-point set contributes \(1+2t\), where \(t\in\{0,1\}\) is the number of chosen points strictly inside the triangle formed by the other three. Summing over all 4-point sets gives

$$Q(m,n)=\binom{N}{4}-C_{\mathrm{bad}}+2\Theta,$$

where \(\Theta\) is the total number of pairs \((T,P)\) such that \(T\) is a nondegenerate lattice triangle and \(P\) is a lattice point strictly inside \(T\).

Step 4: Use Pick's theorem to convert \(\Theta\)

For any lattice triangle, Pick's theorem states

$$A=I+\frac{B}{2}-1,$$

where \(A\) is the area, \(I\) the number of interior lattice points, and \(B\) the number of boundary lattice points. In doubled-area form this becomes

$$2A=2I+B-2.$$

Summing over all nondegenerate lattice triangles yields

$$2\Theta=S_{\mathrm{A}}-S_{\mathrm{B}}+2T,$$

where \(T=\binom{N}{3}-L_3\) is the number of nondegenerate triangles, \(S_{\mathrm{A}}\) is the total doubled area over all nondegenerate triangles, and \(S_{\mathrm{B}}\) is the total boundary-point sum over all nondegenerate triangles.

Step 5: Sum all boundary contributions

For a fixed segment, its contribution to the boundary-point count of every triangle using that segment is exactly \(\gcd(i,j)\). Therefore define

$$S_{\mathrm{seg}}=\sum_{(i,j)\in\mathcal{D}} M_{i,j}g_{i,j}.$$

If we multiply by \(N\), we temporarily pretend that every lattice point can serve as the third vertex. That overcounts points lying on the same maximal lattice line as the segment. If one such line contains \(k\) lattice points, then

$$\sum_{\{A,B\}\subset \ell}\gcd\bigl(|x_B-x_A|,|y_B-y_A|\bigr)=\binom{k+1}{3},$$

because pairs separated by \(r\) primitive steps contribute \(r\), and there are \(k-r\) such pairs. Therefore the forbidden same-line contribution of that line is \(k\binom{k+1}{3}\). Using

$$k\binom{k+1}{3}=2\binom{k}{2}+6\binom{k}{3}+4\binom{k}{4},$$

and summing over all maximal lattice lines, we obtain

$$S_{\mathrm{B}}=N\,S_{\mathrm{seg}}-\left(2\binom{N}{2}+6L_3+4L_4\right).$$

Step 6: Sum all doubled areas by bounding-box size

Every nondegenerate triangle has positive horizontal span and positive vertical span. Group triangles by the size \(i\times j\) of their minimal axis-aligned bounding box, where \(1\le i\le m\) and \(1\le j\le n\). Such a box can be translated in

$$w(i,j)=(m+1-i)(n+1-j)$$

different positions.

For a fixed box size, a direct inclusion-exclusion over all boundary placements of the three vertices simplifies to the exact-box doubled-area sum

$$E(i,j)=\frac{i^2+j^2+11i^2j^2-\gcd(i,j)^2}{3}.$$

This is always an integer. For example, \(E(1,1)=4\) and \(E(2,1)=16\), exactly matching direct enumeration in those boxes. Therefore

$$S_{\mathrm{A}}=\sum_{i=1}^{m}\sum_{j=1}^{n} w(i,j)\,E(i,j).$$

Step 7: Final closed formula and checkpoint

Substituting the previous pieces gives the formula implemented in all three languages:

$$\boxed{Q(m,n)=\binom{N}{4}-\bigl((N-3)L_3-3L_4\bigr)+2\binom{N}{3}-2L_3+S_{\mathrm{A}}-S_{\mathrm{B}}.}$$

For \(m=n=2\), we have \(N=9\), \(L_3=8\), \(L_4=0\), \(S_{\mathrm{A}}=140\), and \(S_{\mathrm{B}}=276\). Hence

$$Q(2,2)=\binom{9}{4}-(6\cdot 8)+2\binom{9}{3}-2\cdot 8+140-276=94,$$

which matches the published checkpoint exactly.

How the Code Works

The C++, Python, and Java implementations evaluate exactly these aggregated sums. They precompute the gcd-derived coefficients, handle horizontal and vertical segments separately, and run the main double loop over positive box widths and heights. The large target is evaluated modulo the prime \(135707531\), while the smaller checkpoints are also verified with exact integer arithmetic.

Complexity Analysis

The dominant work is the double sum over \(1\le i\le m\) and \(1\le j\le n\), so the running time is \(O(mn)\). The auxiliary tables for squares, multiplicities, and gcd-derived coefficients require \(O(m+n)\) memory. Splitting the outer loop across workers improves wall-clock time but does not change the asymptotic complexity.

References

  1. Problem page: https://projecteuler.net/problem=453
  2. Pick's theorem: Wikipedia — Pick's theorem
  3. Lattice points on line segments: Wikipedia — Lattice point
  4. Triangle area by determinant: Wikipedia — Shoelace formula
  5. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle

Problem 453 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

using namespace std;

using int64 = long long;
using i128 = __int128_t;

static constexpr int64 MOD = 135707531;

static int64 mod_add(int64 a, int64 b) {
    a += b;
    if (a >= MOD) a -= MOD;
    return a;
}

static int64 mod_sub(int64 a, int64 b) {
    a -= b;
    if (a < 0) a += MOD;
    return a;
}

static int64 mod_mul(int64 a, int64 b) {
    return static_cast<int64>((static_cast<i128>(a) * b) % MOD);
}

static int64 mod_pow(int64 base, int64 exp) {
    int64 res = 1 % MOD;
    int64 cur = base % MOD;
    while (exp > 0) {
        if (exp & 1) res = mod_mul(res, cur);
        cur = mod_mul(cur, cur);
        exp >>= 1;
    }
    return res;
}

static int64 mod_inv(int64 x) {
    return mod_pow(x, MOD - 2);
}

static string to_string_i128(i128 v) {
    if (v == 0) return "0";
    bool neg = v < 0;
    if (neg) v = -v;
    string s;
    while (v > 0) {
        int digit = static_cast<int>(v % 10);
        s.push_back(static_cast<char>('0' + digit));
        v /= 10;
    }
    if (neg) s.push_back('-');
    reverse(s.begin(), s.end());
    return s;
}

struct ModStats {
    int64 c3 = 0;
    int64 c4 = 0;
    int64 sumg = 0;
    int64 sum_area2 = 0; // Sum of double areas (mod MOD)
};

struct ExactStats {
    i128 c3 = 0;
    i128 c4 = 0;
    i128 sumg = 0;
    i128 sum_area2 = 0; // Sum of double areas (exact)
};

static ModStats compute_stats_mod(int m, int n, int threads) {
    ModStats total;

    const int max_g = max(m, n);
    vector<int64> g_mod(max_g + 1, 0);
    vector<int64> g1_mod(max_g + 1, 0);
    vector<int64> g2_mod(max_g + 1, 0);
    vector<int64> comb2_mod(max_g + 1, 0);
    for (int g = 0; g <= max_g; g++) {
        g_mod[g] = g % MOD;
        g1_mod[g] = (g >= 1) ? (g - 1) % MOD : 0;
        g2_mod[g] = mod_mul(g_mod[g], g_mod[g]);
        if (g >= 2) {
            int64 comb = static_cast<int64>(g - 1) * static_cast<int64>(g - 2) / 2;
            comb2_mod[g] = comb % MOD;
        }
    }

    vector<int64> mx_mod(m + 1, 0);
    vector<int64> ny_mod(n + 1, 0);
    vector<int64> dx2_mod(m + 1, 0);
    vector<int64> dy2_mod(n + 1, 0);

    for (int dx = 0; dx <= m; dx++) {
        mx_mod[dx] = (static_cast<int64>(m + 1 - dx)) % MOD;
        dx2_mod[dx] = mod_mul(dx % MOD, dx % MOD);
    }
    for (int dy = 0; dy <= n; dy++) {
        ny_mod[dy] = (static_cast<int64>(n + 1 - dy)) % MOD;
        dy2_mod[dy] = mod_mul(dy % MOD, dy % MOD);
    }

    const int64 n_plus_1_mod = (n + 1) % MOD;
    const int64 m_plus_1_mod = (m + 1) % MOD;

    // dx > 0, dy = 0
    for (int dx = 1; dx <= m; dx++) {
        int g = dx;
        int64 pairs = mod_mul(n_plus_1_mod, mx_mod[dx]);
        total.sumg = mod_add(total.sumg, mod_mul(pairs, g_mod[g]));
        if (g >= 2) total.c3 = mod_add(total.c3, mod_mul(pairs, g1_mod[g]));
        if (g >= 3) total.c4 = mod_add(total.c4, mod_mul(pairs, comb2_mod[g]));
    }

    // dx = 0, dy > 0
    for (int dy = 1; dy <= n; dy++) {
        int g = dy;
        int64 pairs = mod_mul(m_plus_1_mod, ny_mod[dy]);
        total.sumg = mod_add(total.sumg, mod_mul(pairs, g_mod[g]));
        if (g >= 2) total.c3 = mod_add(total.c3, mod_mul(pairs, g1_mod[g]));
        if (g >= 3) total.c4 = mod_add(total.c4, mod_mul(pairs, comb2_mod[g]));
    }

    int T = threads;
    if (T <= 0) T = static_cast<int>(thread::hardware_concurrency());
    if (T <= 0) T = 4;
    T = min(T, m);
    vector<ModStats> locals(T);
    vector<thread> pool;

    auto worker = [&](int tid, int start_dx, int end_dx) {
        ModStats local;
        for (int dx = start_dx; dx <= end_dx; dx++) {
            int64 mx = mx_mod[dx];
            int64 dx2 = dx2_mod[dx];
            for (int dy = 1; dy <= n; dy++) {
                int g = std::gcd(dx, dy);
                int64 w = mod_mul(mx, ny_mod[dy]);
                int64 pairs = mod_add(w, w); // 2*w

                local.sumg = mod_add(local.sumg, mod_mul(pairs, g_mod[g]));
                if (g >= 2) local.c3 = mod_add(local.c3, mod_mul(pairs, g1_mod[g]));
                if (g >= 3) local.c4 = mod_add(local.c4, mod_mul(pairs, comb2_mod[g]));

                // Sum of double areas uses: (1/3) * sum w*(i^2 + j^2 + 11*i^2*j^2 - gcd^2).
                int64 dy2 = dy2_mod[dy];
                int64 term = mod_add(dx2, dy2);
                term = mod_add(term, mod_mul(11, mod_mul(dx2, dy2)));
                term = mod_sub(term, g2_mod[g]);
                local.sum_area2 = mod_add(local.sum_area2, mod_mul(w, term));
            }
        }
        locals[tid] = local;
    };

    int chunk = (m + T - 1) / T;
    for (int t = 0; t < T; t++) {
        int start_dx = t * chunk + 1;
        int end_dx = min(m, (t + 1) * chunk);
        if (start_dx > end_dx) continue;
        pool.emplace_back(worker, t, start_dx, end_dx);
    }
    for (auto& th : pool) th.join();

    for (const auto& st : locals) {
        total.c3 = mod_add(total.c3, st.c3);
        total.c4 = mod_add(total.c4, st.c4);
        total.sumg = mod_add(total.sumg, st.sumg);
        total.sum_area2 = mod_add(total.sum_area2, st.sum_area2);
    }

    total.sum_area2 = mod_mul(total.sum_area2, mod_inv(3));
    return total;
}

static ExactStats compute_stats_exact(int m, int n) {
    ExactStats total;

    // dx > 0, dy = 0
    for (int dx = 1; dx <= m; dx++) {
        i128 pairs = static_cast<i128>(n + 1) * (m + 1 - dx);
        int g = dx;
        total.sumg += pairs * g;
        if (g >= 2) total.c3 += pairs * (g - 1);
        if (g >= 3) total.c4 += pairs * (g - 1) * (g - 2) / 2;
    }

    // dx = 0, dy > 0
    for (int dy = 1; dy <= n; dy++) {
        i128 pairs = static_cast<i128>(m + 1) * (n + 1 - dy);
        int g = dy;
        total.sumg += pairs * g;
        if (g >= 2) total.c3 += pairs * (g - 1);
        if (g >= 3) total.c4 += pairs * (g - 1) * (g - 2) / 2;
    }

    for (int dx = 1; dx <= m; dx++) {
        for (int dy = 1; dy <= n; dy++) {
            i128 w = static_cast<i128>(m + 1 - dx) * (n + 1 - dy);
            int g = std::gcd(dx, dy);
            i128 pairs = 2 * w;

            total.sumg += pairs * g;
            if (g >= 2) total.c3 += pairs * (g - 1);
            if (g >= 3) total.c4 += pairs * (g - 1) * (g - 2) / 2;

            i128 dx2 = static_cast<i128>(dx) * dx;
            i128 dy2 = static_cast<i128>(dy) * dy;
            i128 term = dx2 + dy2 + 11 * dx2 * dy2 - static_cast<i128>(g) * g;
            total.sum_area2 += w * term;
        }
    }

    total.sum_area2 /= 3;
    return total;
}

static int64 nC2_mod(int64 n) {
    return mod_mul(mod_mul(n, mod_sub(n, 1)), mod_inv(2));
}

static int64 nC3_mod(int64 n) {
    return mod_mul(mod_mul(mod_mul(n, mod_sub(n, 1)), mod_sub(n, 2)), mod_inv(6));
}

static int64 nC4_mod(int64 n) {
    int64 v = mod_mul(n, mod_sub(n, 1));
    v = mod_mul(v, mod_sub(n, 2));
    v = mod_mul(v, mod_sub(n, 3));
    return mod_mul(v, mod_inv(24));
}

static int64 compute_q_mod(int m, int n, int threads) {
    ModStats stats = compute_stats_mod(m, n, threads);

    int64 N = mod_mul((m + 1) % MOD, (n + 1) % MOD);
    int64 C2 = nC2_mod(N);
    int64 C3n = nC3_mod(N);
    int64 C4n = nC4_mod(N);

    int64 S_coll = mod_sub(mod_mul(mod_sub(N, 3), stats.c3), mod_mul(3, stats.c4));

    // LineSum = sum k*C(k+1,3) simplifies to 2*C2 + 6*C3 + 4*C4.
    int64 line_sum = mod_add(mod_mul(2, C2), mod_add(mod_mul(6, stats.c3), mod_mul(4, stats.c4)));

    int64 sumB = mod_sub(mod_mul(N, stats.sumg), line_sum);

    int64 ans = C4n;
    ans = mod_sub(ans, S_coll);
    ans = mod_add(ans, mod_mul(2, C3n));
    ans = mod_sub(ans, mod_mul(2, stats.c3));
    ans = mod_add(ans, stats.sum_area2);
    ans = mod_sub(ans, sumB);
    return ans;
}

static i128 nC2_exact(i128 n) { return n * (n - 1) / 2; }
static i128 nC3_exact(i128 n) { return n * (n - 1) * (n - 2) / 6; }
static i128 nC4_exact(i128 n) { return n * (n - 1) * (n - 2) * (n - 3) / 24; }

static i128 compute_q_exact(int m, int n) {
    ExactStats stats = compute_stats_exact(m, n);
    i128 N = static_cast<i128>(m + 1) * (n + 1);
    i128 C2 = nC2_exact(N);
    i128 C3n = nC3_exact(N);
    i128 C4n = nC4_exact(N);

    i128 S_coll = (N - 3) * stats.c3 - 3 * stats.c4;
    i128 line_sum = 2 * C2 + 6 * stats.c3 + 4 * stats.c4;
    i128 sumB = N * stats.sumg - line_sum;

    i128 ans = C4n - S_coll + 2 * C3n - 2 * stats.c3 + stats.sum_area2 - sumB;
    return ans;
}

static void run_checks() {
    struct Test {
        int m;
        int n;
        i128 expected;
    };
    vector<Test> tests = {
        {2, 2, 94},
        {3, 7, 39590},
        {12, 3, 309000},
        {123, 45, static_cast<i128>(70542215894646LL)}
    };

    cout << "Running validation checks...\n";
    for (const auto& t : tests) {
        i128 got = compute_q_exact(t.m, t.n);
        cout << "Q(" << t.m << ", " << t.n << ") = " << to_string_i128(got)
             << " (Expected: " << to_string_i128(t.expected) << ")";
        if (got == t.expected) {
            cout << " [PASS]\n";
        } else {
            cout << " [FAIL]\n";
        }
        assert(got == t.expected);
    }
    cout << "All validation checks passed!\n";
    cout << "------------------------------------------------\n";
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    run_checks();

    const int m = 12345;
    const int n = 6789;

    int threads = static_cast<int>(thread::hardware_concurrency());
    int64 result = compute_q_mod(m, n, threads);

    cout << "Q(" << m << ", " << n << ") mod " << MOD << " = " << result << "\n";
    return 0;
}

Python

import sys
import math
from multiprocessing import Pool, cpu_count

MOD = 135707531

def mod_pow(base, exp):
    return pow(base, exp, MOD)

def mod_inv(x):
    return mod_pow(x, MOD - 2)

class ModStats:
    def __init__(self):
        self.c3 = 0
        self.c4 = 0
        self.sumg = 0
        self.sum_area2 = 0

def process_chunk(chunk):
    start_dx, end_dx, n, m, g_mod, g1_mod, comb2_mod, g2_mod, mx_mod, ny_mod, dx2_mod, dy2_mod = chunk
    c3 = 0
    c4 = 0
    sumg = 0
    sum_area2 = 0
    
    for dx in range(start_dx, end_dx + 1):
        mx = mx_mod[dx]
        dx2 = dx2_mod[dx]
        for dy in range(1, n + 1):
            g = math.gcd(dx, dy)
            w = (mx * ny_mod[dy]) % MOD
            pairs = (w * 2) % MOD
            
            sumg = (sumg + pairs * g_mod[g]) % MOD
            if g >= 2: c3 = (c3 + pairs * g1_mod[g]) % MOD
            if g >= 3: c4 = (c4 + pairs * comb2_mod[g]) % MOD
            
            dy2 = dy2_mod[dy]
            term = (dx2 + dy2) % MOD
            term = (term + 11 * dx2 * dy2) % MOD
            term = (term - g2_mod[g] + MOD) % MOD
            
            sum_area2 = (sum_area2 + w * term) % MOD
            
    return c3, c4, sumg, sum_area2

def compute_q_mod(m, n):
    max_g = max(m, n)
    g_mod = [g % MOD for g in range(max_g + 1)]
    g1_mod = [max(0, g - 1) % MOD for g in range(max_g + 1)]
    g2_mod = [(g * g) % MOD for g in range(max_g + 1)]
    comb2_mod = [0] * (max_g + 1)
    for g in range(2, max_g + 1):
        comb2_mod[g] = ((g - 1) * (g - 2) // 2) % MOD
        
    mx_mod = [(m + 1 - dx) % MOD for dx in range(m + 1)]
    ny_mod = [(n + 1 - dy) % MOD for dy in range(n + 1)]
    dx2_mod = [(dx % MOD) * (dx % MOD) % MOD for dx in range(m + 1)]
    dy2_mod = [(dy % MOD) * (dy % MOD) % MOD for dy in range(n + 1)]
    
    c3 = 0
    c4 = 0
    sumg = 0
    sum_area2 = 0
    
    n_plus_1_mod = (n + 1) % MOD
    m_plus_1_mod = (m + 1) % MOD
    
    for dx in range(1, m + 1):
        g = dx
        pairs = (n_plus_1_mod * mx_mod[dx]) % MOD
        sumg = (sumg + pairs * g_mod[g]) % MOD
        if g >= 2: c3 = (c3 + pairs * g1_mod[g]) % MOD
        if g >= 3: c4 = (c4 + pairs * comb2_mod[g]) % MOD
        
    for dy in range(1, n + 1):
        g = dy
        pairs = (m_plus_1_mod * ny_mod[dy]) % MOD
        sumg = (sumg + pairs * g_mod[g]) % MOD
        if g >= 2: c3 = (c3 + pairs * g1_mod[g]) % MOD
        if g >= 3: c4 = (c4 + pairs * comb2_mod[g]) % MOD
        
    threads = max(1, cpu_count())
    tasks = []
    chunk_size = max(1, m // threads)
    start_dx = 1
    while start_dx <= m:
        end_dx = min(m, start_dx + chunk_size - 1)
        tasks.append((start_dx, end_dx, n, m, g_mod, g1_mod, comb2_mod, g2_mod, mx_mod, ny_mod, dx2_mod, dy2_mod))
        start_dx = end_dx + 1
        
    with Pool(threads) as pool:
        results = pool.map(process_chunk, tasks)
        
    for res in results:
        c3 = (c3 + res[0]) % MOD
        c4 = (c4 + res[1]) % MOD
        sumg = (sumg + res[2]) % MOD
        sum_area2 = (sum_area2 + res[3]) % MOD
        
    sum_area2 = (sum_area2 * mod_inv(3)) % MOD
    
    N = ((m + 1) % MOD) * ((n + 1) % MOD) % MOD
    C2 = N * (N - 1) % MOD * mod_inv(2) % MOD
    C3n = N * (N - 1) % MOD * (N - 2) % MOD * mod_inv(6) % MOD
    C4n = N * (N - 1) % MOD * (N - 2) % MOD * (N - 3) % MOD * mod_inv(24) % MOD
    
    S_coll = ((N - 3 + MOD) % MOD * c3 - 3 * c4) % MOD
    S_coll = (S_coll + MOD) % MOD
    
    line_sum = (2 * C2 + 6 * c3 + 4 * c4) % MOD
    sumB = (N * sumg - line_sum) % MOD
    sumB = (sumB + MOD) % MOD
    
    ans = C4n
    ans = (ans - S_coll + MOD) % MOD
    ans = (ans + 2 * C3n) % MOD
    ans = (ans - 2 * c3 + MOD) % MOD
    ans = (ans + sum_area2) % MOD
    ans = (ans - sumB + MOD) % MOD
    
    return ans

def solve():
    return str(compute_q_mod(12345, 6789))

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;

public class Euler453 {
    static final long MOD = 135707531;

    static long modPow(long base, long exp) {
        long res = 1;
        long cur = base % MOD;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = (res * cur) % MOD;
            cur = (cur * cur) % MOD;
            exp >>= 1;
        }
        return res;
    }

    static long modInv(long x) {
        return modPow(x, MOD - 2);
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static class ModStats {
        long c3 = 0, c4 = 0, sumg = 0, sumArea2 = 0;
    }

    public static String solve() {
        int m = 12345;
        int n = 6789;

        int maxG = Math.max(m, n);
        long[] gMod = new long[maxG + 1];
        long[] g1Mod = new long[maxG + 1];
        long[] g2Mod = new long[maxG + 1];
        long[] comb2Mod = new long[maxG + 1];

        for (int g = 0; g <= maxG; g++) {
            gMod[g] = g % MOD;
            g1Mod[g] = (g >= 1) ? (g - 1) % MOD : 0;
            g2Mod[g] = (gMod[g] * gMod[g]) % MOD;
            if (g >= 2) {
                long comb = (long) (g - 1) * (g - 2) / 2;
                comb2Mod[g] = comb % MOD;
            }
        }

        long[] mxMod = new long[m + 1];
        long[] nyMod = new long[n + 1];
        long[] dx2Mod = new long[m + 1];
        long[] dy2Mod = new long[n + 1];

        for (int dx = 0; dx <= m; dx++) {
            mxMod[dx] = (m + 1 - dx) % MOD;
            dx2Mod[dx] = ((long) dx * dx) % MOD;
        }
        for (int dy = 0; dy <= n; dy++) {
            nyMod[dy] = (n + 1 - dy) % MOD;
            dy2Mod[dy] = ((long) dy * dy) % MOD;
        }

        ModStats total = new ModStats();
        long nPlus1Mod = (n + 1) % MOD;
        long mPlus1Mod = (m + 1) % MOD;

        for (int dx = 1; dx <= m; dx++) {
            int g = dx;
            long pairs = (nPlus1Mod * mxMod[dx]) % MOD;
            total.sumg = (total.sumg + pairs * gMod[g]) % MOD;
            if (g >= 2)
                total.c3 = (total.c3 + pairs * g1Mod[g]) % MOD;
            if (g >= 3)
                total.c4 = (total.c4 + pairs * comb2Mod[g]) % MOD;
        }

        for (int dy = 1; dy <= n; dy++) {
            int g = dy;
            long pairs = (mPlus1Mod * nyMod[dy]) % MOD;
            total.sumg = (total.sumg + pairs * gMod[g]) % MOD;
            if (g >= 2)
                total.c3 = (total.c3 + pairs * g1Mod[g]) % MOD;
            if (g >= 3)
                total.c4 = (total.c4 + pairs * comb2Mod[g]) % MOD;
        }

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        int chunkSize = m / threads;
        if (chunkSize == 0)
            chunkSize = 1;

        List<int[]> chunks = new ArrayList<>();
        for (int startDx = 1; startDx <= m; startDx += chunkSize) {
            chunks.add(new int[] { startDx, Math.min(m, startDx + chunkSize - 1) });
        }

        long[] sumRes = new long[4];

        chunks.parallelStream().forEach(chunk -> {
            int startDx = chunk[0];
            int endDx = chunk[1];

            long localC3 = 0, localC4 = 0, localSumg = 0, localSumArea2 = 0;

            for (int dx = startDx; dx <= endDx; dx++) {
                long mx = mxMod[dx];
                long dx2 = dx2Mod[dx];
                for (int dy = 1; dy <= n; dy++) {
                    int g = gcd(dx, dy);
                    long w = (mx * nyMod[dy]) % MOD;
                    long pairs = (w * 2) % MOD;

                    localSumg = (localSumg + pairs * gMod[g]) % MOD;
                    if (g >= 2)
                        localC3 = (localC3 + pairs * g1Mod[g]) % MOD;
                    if (g >= 3)
                        localC4 = (localC4 + pairs * comb2Mod[g]) % MOD;

                    long dy2 = dy2Mod[dy];
                    long term = (dx2 + dy2) % MOD;
                    term = (term + 11 * dx2 % MOD * dy2) % MOD;
                    term = (term - g2Mod[g] + MOD) % MOD;

                    localSumArea2 = (localSumArea2 + w * term) % MOD;
                }
            }

            synchronized (sumRes) {
                sumRes[0] = (sumRes[0] + localC3) % MOD;
                sumRes[1] = (sumRes[1] + localC4) % MOD;
                sumRes[2] = (sumRes[2] + localSumg) % MOD;
                sumRes[3] = (sumRes[3] + localSumArea2) % MOD;
            }
        });

        total.c3 = (total.c3 + sumRes[0]) % MOD;
        total.c4 = (total.c4 + sumRes[1]) % MOD;
        total.sumg = (total.sumg + sumRes[2]) % MOD;
        total.sumArea2 = (total.sumArea2 + sumRes[3]) % MOD;

        total.sumArea2 = (total.sumArea2 * modInv(3)) % MOD;

        long bigN = (mPlus1Mod * nPlus1Mod) % MOD;
        long C2 = bigN * (bigN - 1) % MOD * modInv(2) % MOD;
        long C3n = bigN * (bigN - 1) % MOD * (bigN - 2 + MOD) % MOD * modInv(6) % MOD;
        long C4n = bigN * (bigN - 1) % MOD * (bigN - 2 + MOD) % MOD * (bigN - 3 + MOD) % MOD * modInv(24) % MOD;

        long SColl = ((bigN - 3 + MOD) % MOD * total.c3 % MOD - 3 * total.c4 % MOD + MOD) % MOD;

        long lineSum = (2 * C2 + 6 * total.c3 + 4 * total.c4) % MOD;
        long sumB = (bigN * total.sumg % MOD - lineSum + MOD) % MOD;

        long ans = C4n;
        ans = (ans - SColl + MOD) % MOD;
        ans = (ans + 2 * C3n) % MOD;
        ans = (ans - 2 * total.c3 % MOD + MOD) % MOD;
        ans = (ans + total.sumArea2) % MOD;
        ans = (ans - sumB + MOD) % MOD;

        return Long.toString(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}