Problem 994: Counting Triangles

Project Euler Problem 994 Solution

EulerSolve provides an optimized solution for Project Euler Problem 994, Counting Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every pair $1\le i\le m$, $1\le j\le n$, the picture contains the segment from the lower point $(i,1)$ to the upper point $(j,2)$. A triangle is formed by three of these segments when the three pairwise intersections are real and not all the same point. Intersections at shared endpoints count, and intersections created inside other segments count as well. The target is $T(1234\cdot 10^8,2345\cdot 10^8)$ modulo $10^9+7$. The given checks are $T(2,3)=8$, $T(3,5)=146$, and $T(12,23)=756716$. A geometric drawing is far too large, so the implementation counts triples of segments algebraically. Mathematical Approach 1. Dual view of a segment Represent the segment from $(i,1)$ to $(j,2)$ by the grid point $(i,j)$ in an $m\times n$ rectangle. Two distinct segments $(i_1,j_1)$ and $(i_2,j_2)$ intersect exactly when $$i_1=i_2,\qquad\text{or}\qquad j_1=j_2,\qquad\text{or}\qquad (i_1-i_2)(j_1-j_2)\lt 0.$$ The first two cases are shared lower or upper endpoints. The third case says that the two endpoints appear in opposite order on the two horizontal lines, so the segments cross in the interior. 2. Counting triples that intersect pairwise We first count triples of segments for which every pair intersects. Some of these will later be removed because all three meet at one point....

Detailed mathematical approach

Problem Summary

For every pair $1\le i\le m$, $1\le j\le n$, the picture contains the segment from the lower point $(i,1)$ to the upper point $(j,2)$. A triangle is formed by three of these segments when the three pairwise intersections are real and not all the same point. Intersections at shared endpoints count, and intersections created inside other segments count as well.

The target is $T(1234\cdot 10^8,2345\cdot 10^8)$ modulo $10^9+7$. The given checks are $T(2,3)=8$, $T(3,5)=146$, and $T(12,23)=756716$. A geometric drawing is far too large, so the implementation counts triples of segments algebraically.

Mathematical Approach

1. Dual view of a segment

Represent the segment from $(i,1)$ to $(j,2)$ by the grid point $(i,j)$ in an $m\times n$ rectangle. Two distinct segments $(i_1,j_1)$ and $(i_2,j_2)$ intersect exactly when

$$i_1=i_2,\qquad\text{or}\qquad j_1=j_2,\qquad\text{or}\qquad (i_1-i_2)(j_1-j_2)\lt 0.$$

The first two cases are shared lower or upper endpoints. The third case says that the two endpoints appear in opposite order on the two horizontal lines, so the segments cross in the interior.

2. Counting triples that intersect pairwise

We first count triples of segments for which every pair intersects. Some of these will later be removed because all three meet at one point. The count splits cleanly by how many lower endpoints and upper endpoints are used.

If the triple uses three different lower endpoints and three different upper endpoints, then after the lower endpoints are ordered increasingly, the upper endpoints must be ordered decreasingly. Thus there is one valid matching for every choice of three lower and three upper endpoints:

$$\binom m3\binom n3.$$

If it uses three lower endpoints but only two upper endpoints, choose the three lower endpoints and an ordered pair of distinct upper endpoints; the singleton upper endpoint must be attached to the extreme lower endpoint on the side that makes it cross the other two. This gives $\binom m3 n(n-1)$. By symmetry, the case of two lower endpoints and three upper endpoints contributes $m(m-1)\binom n3$.

Finally, in a $2\times 2$ rectangle of endpoints, exactly two of the four three-corner choices are pairwise intersecting. Therefore the $2$-by-$2$ case contributes

$$2\binom m2\binom n2={m(m-1)n(n-1)\over 2}.$$

So the pairwise-intersecting candidate count is

$$P(m,n)=\binom m3\binom n3+\binom m3 n(n-1)+m(m-1)\binom n3+{m(m-1)n(n-1)\over 2}.$$

3. Removing triples concurrent at one point

A candidate triple fails to be a triangle precisely when the three segment-lines are concurrent. In the dual grid this is the same as the three points $(i,j)$ being collinear on a line of negative slope.

Take the two extreme dual points of such a collinear triple. Let their horizontal separation be $A$ and their vertical separation be $B$, with both positive. There are $(m-A)(n-B)$ ways to place those two extremes in the rectangle. The number of possible middle lattice points on the open segment between them is

$$\gcd(A,B)-1.$$

Hence the number of concurrent triples is

$$C(m,n)=\sum_{A=1}^{m-1}\sum_{B=1}^{n-1}(\gcd(A,B)-1)(m-A)(n-B).$$

4. Totient transform

The double sum above is still too large. Use the classical identity

$$\gcd(A,B)=\sum_{d\mid\gcd(A,B)}\varphi(d),$$

so that

$$\gcd(A,B)-1=\sum_{\substack{d\mid A,\ d\mid B\\d\ge 2}}\varphi(d).$$

Switching the order of summation gives

$$C(m,n)=\sum_{d=2}^{\min(m,n)-1}\varphi(d) \left(\sum_{r=1}^{\lfloor(m-1)/d\rfloor}(m-dr)\right) \left(\sum_{s=1}^{\lfloor(n-1)/d\rfloor}(n-ds)\right).$$

Define

$$Q_m(d)=\left\lfloor {m-1\over d}\right\rfloor,\qquad A_m(d)=mQ_m(d)-d{Q_m(d)(Q_m(d)+1)\over 2}.$$

Then $C(m,n)=\sum_d \varphi(d)A_m(d)A_n(d)$.

5. Summatory totients and quotient blocks

The values $Q_m(d)$ and $Q_n(d)$ are constant on long intervals of $d$. On such an interval, $A_m(d)A_n(d)$ is a quadratic polynomial in $d$, so only three prefix sums are required:

$$\Phi_k(x)=\sum_{d\le x}d^k\varphi(d),\qquad k=0,1,2.$$

The code builds these sums up to $5\cdot 10^6$ by a linear sieve and evaluates larger arguments with the divisor-summatory recurrence

$$\Phi_k(n)=\sum_{r=1}^{n}r^{k+1} -\sum_{\ell=2}^{n}\left(\sum_{d=\ell}^{h}d^k\right)\Phi_k\!\left(\left\lfloor {n\over \ell}\right\rfloor\right),$$

where $h=\left\lfloor n/\left\lfloor n/\ell\right\rfloor\right\rfloor$. Grouping equal quotients makes this recurrence fast enough for the huge target.

6. Final formula and checks

The answer is

$$\boxed{T(m,n)=P(m,n)-C(m,n)\pmod {10^9+7}}.$$

For $m=2,n=3$, the concurrent sum is zero and $P(2,3)=8$, so $T(2,3)=8$. For $m=3,n=5$, $P(3,5)=150$. The only nonzero concurrent block is $d=2$, which contributes $4$, giving $150-4=146$. The larger checkpoint $T(12,23)=756716$ is also asserted in the implementation.

How the Code Works

The C++ program first initializes modular arithmetic helpers, binomial forms, and power-sum formulas up to degree three. pairwise_triangle_candidates evaluates $P(m,n)$ directly.

TotientSummatory builds $\varphi(d)$, $\Phi_0$, $\Phi_1$, and $\Phi_2$ up to the sieve limit, then memoizes the recursive prefix calls for larger $n$. concurrent_triples walks over quotient blocks of $d$, expands $A_m(d)A_n(d)$, and subtracts the required totient-weighted sum. The small brute-force checker enumerates all segment triples for $1\le m,n\le 5$, then verifies the three official checkpoints.

Complexity Analysis

The fixed sieve costs $O(L)$ time and memory for $L=5\cdot 10^6$. The main summation uses quotient blocks, so the number of outer iterations is $O(\sqrt{\min(m,n)})$, with memoized summatory-totient queries. The memory use is dominated by the three prefix arrays and the totient table.

A direct enumeration of all triples of the $mn$ segments would be impossible. The formula replaces that $O((mn)^3)$ geometry with arithmetic over divisors and floor-division blocks.

Footnotes and References

Problem page: Project Euler 994
Euler's totient function: Wikipedia - Euler's totient function
Floor-sum decomposition: Wikipedia - Dirichlet hyperbola method
Line arrangement: Wikipedia - Arrangement of lines

Problem 994 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i64 = std::int64_t;
using i128 = __int128_t;

constexpr u64 TARGET_M = 1234ULL * 100'000'000ULL;
constexpr u64 TARGET_N = 2345ULL * 100'000'000ULL;
constexpr u32 MOD = 1'000'000'007U;
constexpr int SIEVE_LIMIT = 5'000'000;
constexpr u64 INV2 = 500'000'004ULL;
constexpr u64 INV3 = 333'333'336ULL;
constexpr u64 INV6 = 166'666'668ULL;

u32 add_mod(const u32 a, const u32 b) {
    const u32 s = a + b;
    return s >= MOD ? s - MOD : s;
}

u32 sub_mod(const u32 a, const u32 b) {
    return a >= b ? a - b : static_cast<u32>(static_cast<u64>(a) + MOD - b);
}

u32 mul_mod(const u64 a, const u64 b) {
    return static_cast<u32>((a % MOD) * (b % MOD) % MOD);
}

u32 sum_power(const int power, const u64 n) {
    const u64 a = n % MOD;
    const u64 b = (n + 1) % MOD;
    if (power == 0) {
        return static_cast<u32>(a);
    }
    if (power == 1) {
        return static_cast<u32>(a * b % MOD * INV2 % MOD);
    }
    if (power == 2) {
        const u64 c = (2 * (n % MOD) + 1) % MOD;
        return static_cast<u32>(a * b % MOD * c % MOD * INV6 % MOD);
    }
    const u64 t = a * b % MOD * INV2 % MOD;
    return static_cast<u32>(t * t % MOD);
}

u32 range_power(const int power, const u64 lo, const u64 hi) {
    return sub_mod(sum_power(power, hi), sum_power(power, lo - 1));
}

u32 choose2(const u64 n) {
    const u64 a = n % MOD;
    const u64 b = (n - 1) % MOD;
    return static_cast<u32>(a * b % MOD * INV2 % MOD);
}

u32 choose3(const u64 n) {
    if (n < 3) {
        return 0;
    }
    const u64 a = n % MOD;
    const u64 b = (n - 1) % MOD;
    const u64 c = (n - 2) % MOD;
    return static_cast<u32>(a * b % MOD * c % MOD * INV6 % MOD);
}

u32 triangle_number(const u64 n) {
    const u64 a = n % MOD;
    const u64 b = (n + 1) % MOD;
    return static_cast<u32>(a * b % MOD * INV2 % MOD);
}

struct TotientSummatory {
    explicit TotientSummatory(const int limit) : limit(limit) {
        build();
        for (auto& table : memo) {
            table.reserve(200'000);
        }
    }

    u32 prefix(const int power, const u64 n) {
        if (n <= static_cast<u64>(limit)) {
            return sums[static_cast<std::size_t>(power)][static_cast<std::size_t>(n)];
        }

        auto& table = memo[static_cast<std::size_t>(power)];
        const auto it = table.find(n);
        if (it != table.end()) {
            return it->second;
        }

        u32 result = sum_power(power + 1, n);
        for (u64 lo = 2; lo <= n;) {
            const u64 q = n / lo;
            const u64 hi = n / q;
            const u32 block = range_power(power, lo, hi);
            result = sub_mod(result, mul_mod(block, prefix(power, q)));
            lo = hi + 1;
        }

        table.emplace(n, result);
        return result;
    }

private:
    int limit;
    std::array<std::vector<u32>, 3> sums;
    std::array<std::unordered_map<u64, u32>, 3> memo;

    void build() {
        std::vector<int> phi(static_cast<std::size_t>(limit + 1), 0);
        std::vector<int> primes;
        std::vector<unsigned char> composite(static_cast<std::size_t>(limit + 1), 0);
        phi[1] = 1;

        for (int i = 2; i <= limit; ++i) {
            if (!composite[static_cast<std::size_t>(i)]) {
                primes.push_back(i);
                phi[static_cast<std::size_t>(i)] = i - 1;
            }
            for (const int p : primes) {
                const i64 v = static_cast<i64>(i) * p;
                if (v > limit) {
                    break;
                }
                composite[static_cast<std::size_t>(v)] = 1;
                if (i % p == 0) {
                    phi[static_cast<std::size_t>(v)] = phi[static_cast<std::size_t>(i)] * p;
                    break;
                }
                phi[static_cast<std::size_t>(v)] = phi[static_cast<std::size_t>(i)] * (p - 1);
            }
        }

        for (auto& row : sums) {
            row.assign(static_cast<std::size_t>(limit + 1), 0);
        }
        for (int i = 1; i <= limit; ++i) {
            const u64 im = static_cast<u64>(i) % MOD;
            const u64 ph = static_cast<u64>(phi[static_cast<std::size_t>(i)]) % MOD;
            sums[0][static_cast<std::size_t>(i)] = add_mod(sums[0][static_cast<std::size_t>(i - 1)],
                                                           static_cast<u32>(ph));
            sums[1][static_cast<std::size_t>(i)] = add_mod(sums[1][static_cast<std::size_t>(i - 1)],
                                                           static_cast<u32>(im * ph % MOD));
            sums[2][static_cast<std::size_t>(i)] = add_mod(sums[2][static_cast<std::size_t>(i - 1)],
                                                           static_cast<u32>(im * im % MOD * ph % MOD));
        }
    }
};

u32 pairwise_triangle_candidates(const u64 m, const u64 n) {
    u32 total = mul_mod(choose3(m), choose3(n));
    total = add_mod(total, mul_mod(mul_mod(choose3(m), n % MOD), (n - 1) % MOD));
    total = add_mod(total, mul_mod(mul_mod(m % MOD, (m - 1) % MOD), choose3(n)));
    const u32 corner = mul_mod(mul_mod(m % MOD, (m - 1) % MOD), mul_mod(n % MOD, (n - 1) % MOD));
    total = add_mod(total, mul_mod(corner, INV2));
    return total;
}

u32 concurrent_triples(const u64 m, const u64 n, TotientSummatory& totient) {
    const u64 limit = std::min(m, n) - 1;
    u32 total = 0;

    for (u64 lo = 2; lo <= limit;) {
        const u64 qm = (m - 1) / lo;
        const u64 qn = (n - 1) / lo;
        const u64 hi = std::min((m - 1) / qm, (n - 1) / qn);

        const u32 s0 = sub_mod(totient.prefix(0, hi), totient.prefix(0, lo - 1));
        const u32 s1 = sub_mod(totient.prefix(1, hi), totient.prefix(1, lo - 1));
        const u32 s2 = sub_mod(totient.prefix(2, hi), totient.prefix(2, lo - 1));

        const u32 am0 = mul_mod(qm % MOD, m % MOD);
        const u32 an0 = mul_mod(qn % MOD, n % MOD);
        const u32 am1 = triangle_number(qm);
        const u32 an1 = triangle_number(qn);

        u32 term = mul_mod(mul_mod(am0, an0), s0);
        term = sub_mod(term, mul_mod(add_mod(mul_mod(am0, an1), mul_mod(am1, an0)), s1));
        term = add_mod(term, mul_mod(mul_mod(am1, an1), s2));
        total = add_mod(total, term);

        lo = hi + 1;
    }

    return total;
}

u32 solve_mod(const u64 m, const u64 n, TotientSummatory& totient) {
    return sub_mod(pairwise_triangle_candidates(m, n), concurrent_triples(m, n, totient));
}

bool pair_intersects(const std::pair<int, int>& a, const std::pair<int, int>& b) {
    return a.first == b.first || a.second == b.second ||
           static_cast<i64>(a.first - b.first) * static_cast<i64>(a.second - b.second) < 0;
}

bool collinear(const std::pair<int, int>& a, const std::pair<int, int>& b, const std::pair<int, int>& c) {
    return static_cast<i64>(b.first - a.first) * static_cast<i64>(c.second - a.second) ==
           static_cast<i64>(b.second - a.second) * static_cast<i64>(c.first - a.first);
}

u64 brute_count(const int m, const int n) {
    std::vector<std::pair<int, int>> segments;
    for (int i = 1; i <= m; ++i) {
        for (int j = 1; j <= n; ++j) {
            segments.emplace_back(i, j);
        }
    }

    u64 total = 0;
    for (std::size_t a = 0; a < segments.size(); ++a) {
        for (std::size_t b = a + 1; b < segments.size(); ++b) {
            for (std::size_t c = b + 1; c < segments.size(); ++c) {
                if (pair_intersects(segments[a], segments[b]) &&
                    pair_intersects(segments[a], segments[c]) &&
                    pair_intersects(segments[b], segments[c]) &&
                    !collinear(segments[a], segments[b], segments[c])) {
                    ++total;
                }
            }
        }
    }
    return total;
}

void run_checkpoints(TotientSummatory& totient) {
    for (int m = 1; m <= 5; ++m) {
        for (int n = 1; n <= 5; ++n) {
            assert(solve_mod(m, n, totient) == brute_count(m, n));
        }
    }

    assert(solve_mod(2, 3, totient) == 8U);
    assert(solve_mod(3, 5, totient) == 146U);
    assert(solve_mod(12, 23, totient) == 756'716U);
}

}  // namespace

int main(int argc, char** argv) {
    bool should_run_checkpoints = true;
    for (int i = 1; i < argc; ++i) {
        if (std::string(argv[i]) == "--skip-checkpoints") {
            should_run_checkpoints = false;
        }
    }

    TotientSummatory totient(SIEVE_LIMIT);
    if (should_run_checkpoints) {
        run_checkpoints(totient);
    }

    std::cout << solve_mod(TARGET_M, TARGET_N, totient) << '\n';
    return 0;
}

Python

import sys
from array import array

TARGET_M = 1234 * 100_000_000
TARGET_N = 2345 * 100_000_000
MOD = 1_000_000_007
SIEVE_LIMIT = 5_000_000
INV2 = 500_000_004
INV3 = 333_333_336
INV6 = 166_666_668


def add_mod(a, b):
    s = a + b
    return s - MOD if s >= MOD else s


def sub_mod(a, b):
    return a - b if a >= b else a + MOD - b


def mul_mod(a, b):
    return (a % MOD) * (b % MOD) % MOD


def sum_power(power, n):
    a = n % MOD
    b = (n + 1) % MOD
    if power == 0:
        return a
    if power == 1:
        return a * b % MOD * INV2 % MOD
    if power == 2:
        c = (2 * (n % MOD) + 1) % MOD
        return a * b % MOD * c % MOD * INV6 % MOD
    t = a * b % MOD * INV2 % MOD
    return t * t % MOD


def range_power(power, lo, hi):
    return sub_mod(sum_power(power, hi), sum_power(power, lo - 1))


def choose2(n):
    return n % MOD * ((n - 1) % MOD) % MOD * INV2 % MOD


def choose3(n):
    if n < 3:
        return 0
    return n % MOD * ((n - 1) % MOD) % MOD * ((n - 2) % MOD) % MOD * INV6 % MOD


def triangle_number(n):
    return n % MOD * ((n + 1) % MOD) % MOD * INV2 % MOD


class TotientSummatory:
    def __init__(self, limit):
        self.limit = limit
        self.sums = [array("I", [0]) * (limit + 1) for _ in range(3)]
        self.memo = [{} for _ in range(3)]
        self._build()

    def prefix(self, power, n):
        if n <= self.limit:
            return self.sums[power][n]

        table = self.memo[power]
        cached = table.get(n)
        if cached is not None:
            return cached

        result = sum_power(power + 1, n)
        lo = 2
        while lo <= n:
            q = n // lo
            hi = n // q
            block = range_power(power, lo, hi)
            result = (result - block * self.prefix(power, q)) % MOD
            lo = hi + 1

        table[n] = result
        return result

    def _build(self):
        limit = self.limit
        phi = array("I", range(limit + 1))
        phi[1] = 1

        for i in range(2, limit + 1):
            if phi[i] == i:
                for j in range(i, limit + 1, i):
                    phi[j] -= phi[j] // i

        s0, s1, s2 = self.sums
        for i in range(1, limit + 1):
            im = i % MOD
            ph = phi[i] % MOD
            s0[i] = (s0[i - 1] + ph) % MOD
            s1[i] = (s1[i - 1] + im * ph) % MOD
            s2[i] = (s2[i - 1] + im * im % MOD * ph) % MOD


def pairwise_triangle_candidates(m, n):
    total = mul_mod(choose3(m), choose3(n))
    total = add_mod(total, mul_mod(mul_mod(choose3(m), n % MOD), (n - 1) % MOD))
    total = add_mod(total, mul_mod(mul_mod(m % MOD, (m - 1) % MOD), choose3(n)))
    corner = mul_mod(mul_mod(m % MOD, (m - 1) % MOD), mul_mod(n % MOD, (n - 1) % MOD))
    total = add_mod(total, mul_mod(corner, INV2))
    return total


def concurrent_triples(m, n, totient):
    limit = min(m, n) - 1
    total = 0
    lo = 2

    while lo <= limit:
        qm = (m - 1) // lo
        qn = (n - 1) // lo
        hi = min((m - 1) // qm, (n - 1) // qn)

        s0 = sub_mod(totient.prefix(0, hi), totient.prefix(0, lo - 1))
        s1 = sub_mod(totient.prefix(1, hi), totient.prefix(1, lo - 1))
        s2 = sub_mod(totient.prefix(2, hi), totient.prefix(2, lo - 1))

        am0 = mul_mod(qm, m)
        an0 = mul_mod(qn, n)
        am1 = triangle_number(qm)
        an1 = triangle_number(qn)

        term = mul_mod(mul_mod(am0, an0), s0)
        term = sub_mod(term, mul_mod(add_mod(mul_mod(am0, an1), mul_mod(am1, an0)), s1))
        term = add_mod(term, mul_mod(mul_mod(am1, an1), s2))
        total = add_mod(total, term)

        lo = hi + 1

    return total


def solve_mod(m, n, totient):
    return sub_mod(pairwise_triangle_candidates(m, n), concurrent_triples(m, n, totient))


def pair_intersects(a, b):
    return a[0] == b[0] or a[1] == b[1] or (a[0] - b[0]) * (a[1] - b[1]) < 0


def collinear(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) == (b[1] - a[1]) * (c[0] - a[0])


def brute_count(m, n):
    segments = [(i, j) for i in range(1, m + 1) for j in range(1, n + 1)]
    total = 0

    for a in range(len(segments)):
        for b in range(a + 1, len(segments)):
            for c in range(b + 1, len(segments)):
                sa = segments[a]
                sb = segments[b]
                sc = segments[c]
                if (
                    pair_intersects(sa, sb)
                    and pair_intersects(sa, sc)
                    and pair_intersects(sb, sc)
                    and not collinear(sa, sb, sc)
                ):
                    total += 1

    return total


def run_checkpoints(totient):
    for m in range(1, 6):
        for n in range(1, 6):
            assert solve_mod(m, n, totient) == brute_count(m, n)

    assert solve_mod(2, 3, totient) == 8
    assert solve_mod(3, 5, totient) == 146
    assert solve_mod(12, 23, totient) == 756_716


def main(argv=None):
    argv = sys.argv if argv is None else argv
    should_run_checkpoints = True

    for arg in argv[1:]:
        if arg == "--skip-checkpoints":
            should_run_checkpoints = False
            continue
        print(f"Unknown argument: {arg}", file=sys.stderr)
        return 1

    totient = TotientSummatory(SIEVE_LIMIT)
    if should_run_checkpoints:
        run_checkpoints(totient)

    print(solve_mod(TARGET_M, TARGET_N, totient))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())

Java

import java.util.HashMap;

public class Euler994 {
    private static final long TARGET_M = 1234L * 100_000_000L;
    private static final long TARGET_N = 2345L * 100_000_000L;
    private static final int MOD = 1_000_000_007;
    private static final int SIEVE_LIMIT = 5_000_000;
    private static final long INV2 = 500_000_004L;
    private static final long INV3 = 333_333_336L;
    private static final long INV6 = 166_666_668L;

    private static int addMod(int a, int b) {
        int s = a + b;
        return s >= MOD ? s - MOD : s;
    }

    private static int subMod(int a, int b) {
        return a >= b ? a - b : (int) ((long) a + MOD - b);
    }

    private static int mulMod(long a, long b) {
        return (int) ((a % MOD) * (b % MOD) % MOD);
    }

    private static int sumPower(int power, long n) {
        long a = n % MOD;
        long b = (n + 1) % MOD;
        if (power == 0) {
            return (int) a;
        }
        if (power == 1) {
            return (int) (a * b % MOD * INV2 % MOD);
        }
        if (power == 2) {
            long c = (2 * (n % MOD) + 1) % MOD;
            return (int) (a * b % MOD * c % MOD * INV6 % MOD);
        }
        long t = a * b % MOD * INV2 % MOD;
        return (int) (t * t % MOD);
    }

    private static int rangePower(int power, long lo, long hi) {
        return subMod(sumPower(power, hi), sumPower(power, lo - 1));
    }

    private static int choose2(long n) {
        long a = n % MOD;
        long b = (n - 1) % MOD;
        return (int) (a * b % MOD * INV2 % MOD);
    }

    private static int choose3(long n) {
        if (n < 3L) {
            return 0;
        }
        long a = n % MOD;
        long b = (n - 1) % MOD;
        long c = (n - 2) % MOD;
        return (int) (a * b % MOD * c % MOD * INV6 % MOD);
    }

    private static int triangleNumber(long n) {
        long a = n % MOD;
        long b = (n + 1) % MOD;
        return (int) (a * b % MOD * INV2 % MOD);
    }

    private static final class TotientSummatory {
        private final int limit;
        private final int[][] sums = new int[3][];
        private final HashMap<Long, Integer>[] memo;

        @SuppressWarnings("unchecked")
        private TotientSummatory(int limit) {
            this.limit = limit;
            for (int i = 0; i < 3; ++i) {
                sums[i] = new int[limit + 1];
            }
            memo = new HashMap[3];
            for (int i = 0; i < 3; ++i) {
                memo[i] = new HashMap<>(262_144);
            }
            build();
        }

        private int prefix(int power, long n) {
            if (n <= limit) {
                return sums[power][(int) n];
            }

            Integer cached = memo[power].get(n);
            if (cached != null) {
                return cached;
            }

            int result = sumPower(power + 1, n);
            for (long lo = 2L; lo <= n;) {
                long q = n / lo;
                long hi = n / q;
                int block = rangePower(power, lo, hi);
                result = subMod(result, mulMod(block, prefix(power, q)));
                lo = hi + 1;
            }

            memo[power].put(n, result);
            return result;
        }

        private void build() {
            int[] phi = new int[limit + 1];
            int[] primes = new int[limit];
            int primeCount = 0;
            boolean[] composite = new boolean[limit + 1];
            phi[1] = 1;

            for (int i = 2; i <= limit; ++i) {
                if (!composite[i]) {
                    primes[primeCount++] = i;
                    phi[i] = i - 1;
                }

                for (int j = 0; j < primeCount; ++j) {
                    int p = primes[j];
                    long v = (long) i * p;
                    if (v > limit) {
                        break;
                    }
                    composite[(int) v] = true;
                    if (i % p == 0) {
                        phi[(int) v] = phi[i] * p;
                        break;
                    }
                    phi[(int) v] = phi[i] * (p - 1);
                }
            }

            for (int i = 1; i <= limit; ++i) {
                long im = i % (long) MOD;
                long ph = phi[i] % (long) MOD;
                sums[0][i] = addMod(sums[0][i - 1], (int) ph);
                sums[1][i] = addMod(sums[1][i - 1], (int) (im * ph % MOD));
                sums[2][i] = addMod(sums[2][i - 1], (int) (im * im % MOD * ph % MOD));
            }
        }
    }

    private static int pairwiseTriangleCandidates(long m, long n) {
        int total = mulMod(choose3(m), choose3(n));
        total = addMod(total, mulMod(mulMod(choose3(m), n % MOD), (n - 1) % MOD));
        total = addMod(total, mulMod(mulMod(m % MOD, (m - 1) % MOD), choose3(n)));
        int corner = mulMod(mulMod(m % MOD, (m - 1) % MOD), mulMod(n % MOD, (n - 1) % MOD));
        total = addMod(total, mulMod(corner, INV2));
        return total;
    }

    private static int concurrentTriples(long m, long n, TotientSummatory totient) {
        long limit = Math.min(m, n) - 1L;
        int total = 0;

        for (long lo = 2L; lo <= limit;) {
            long qm = (m - 1) / lo;
            long qn = (n - 1) / lo;
            long hi = Math.min((m - 1) / qm, (n - 1) / qn);

            int s0 = subMod(totient.prefix(0, hi), totient.prefix(0, lo - 1));
            int s1 = subMod(totient.prefix(1, hi), totient.prefix(1, lo - 1));
            int s2 = subMod(totient.prefix(2, hi), totient.prefix(2, lo - 1));

            int am0 = mulMod(qm, m);
            int an0 = mulMod(qn, n);
            int am1 = triangleNumber(qm);
            int an1 = triangleNumber(qn);

            int term = mulMod(mulMod(am0, an0), s0);
            term = subMod(term, mulMod(addMod(mulMod(am0, an1), mulMod(am1, an0)), s1));
            term = addMod(term, mulMod(mulMod(am1, an1), s2));
            total = addMod(total, term);

            lo = hi + 1;
        }

        return total;
    }

    private static int solveMod(long m, long n, TotientSummatory totient) {
        return subMod(pairwiseTriangleCandidates(m, n), concurrentTriples(m, n, totient));
    }

    private static boolean pairIntersects(int[] a, int[] b) {
        return a[0] == b[0] || a[1] == b[1] || (long) (a[0] - b[0]) * (a[1] - b[1]) < 0L;
    }

    private static boolean collinear(int[] a, int[] b, int[] c) {
        return (long) (b[0] - a[0]) * (c[1] - a[1])
            == (long) (b[1] - a[1]) * (c[0] - a[0]);
    }

    private static long bruteCount(int m, int n) {
        int[][] segments = new int[m * n][2];
        int size = 0;
        for (int i = 1; i <= m; ++i) {
            for (int j = 1; j <= n; ++j) {
                segments[size][0] = i;
                segments[size][1] = j;
                ++size;
            }
        }

        long total = 0L;
        for (int a = 0; a < size; ++a) {
            for (int b = a + 1; b < size; ++b) {
                for (int c = b + 1; c < size; ++c) {
                    if (pairIntersects(segments[a], segments[b])
                        && pairIntersects(segments[a], segments[c])
                        && pairIntersects(segments[b], segments[c])
                        && !collinear(segments[a], segments[b], segments[c])) {
                        ++total;
                    }
                }
            }
        }
        return total;
    }

    private static void require(boolean condition, String message) {
        if (!condition) {
            throw new IllegalStateException(message);
        }
    }

    private static void runCheckpoints(TotientSummatory totient) {
        for (int m = 1; m <= 5; ++m) {
            for (int n = 1; n <= 5; ++n) {
                require(solveMod(m, n, totient) == bruteCount(m, n), "Brute checkpoint failed");
            }
        }

        require(solveMod(2L, 3L, totient) == 8, "Checkpoint failed for T(2,3)");
        require(solveMod(3L, 5L, totient) == 146, "Checkpoint failed for T(3,5)");
        require(solveMod(12L, 23L, totient) == 756_716, "Checkpoint failed for T(12,23)");
    }

    public static void main(String[] args) {
        boolean shouldRunCheckpoints = true;
        for (String arg : args) {
            if ("--skip-checkpoints".equals(arg)) {
                shouldRunCheckpoints = false;
                continue;
            }
            System.err.println("Unknown argument: " + arg);
            System.exit(1);
        }

        TotientSummatory totient = new TotientSummatory(SIEVE_LIMIT);
        if (shouldRunCheckpoints) {
            runCheckpoints(totient);
        }

        System.out.println(solveMod(TARGET_M, TARGET_N, totient));
    }
}