Problem 883: Remarkable Triangles

Project Euler Problem 883 Solution

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

Problem Summary Problem 883 asks for the value $T(R_2)$ for remarkable triangles at $R_2=2{,}000{,}000$. The implementations reduce the geometry to counting points in the triangular lattice, where a vector written in basis coordinates $(u,v)$ has squared length $$Q(u,v)=u^2+uv+v^2.$$ The final answer is not obtained by one direct geometric count: it is a weighted sum over many scaled lattice counts, with Möbius inversion removing non-primitive overcounting. Mathematical Approach The triangular-lattice geometry behind remarkable triangles is encoded by the quadratic form $Q(u,v)=u^2+uv+v^2$. The problem therefore becomes an arithmetic-geometric counting problem built from this norm form. Step 1: Convert the Geometry into Lattice Counts For a threshold $M$, define $$A(M)=\#\{(u,v)\in \mathbb{Z}^2:Q(u,v)\le M\},$$ $$A_{\equiv}(M)=\#\{(u,v)\in \mathbb{Z}^2:Q(u,v)\le M,\ u\equiv v\pmod{3}\}.$$ The unrestricted count $A(M)$ measures all admissible lattice vectors up to squared length $M$. The restricted count $A_{\equiv}(M)$ keeps only the congruence class needed when a factor of $3$ has to be absorbed separately....

Detailed mathematical approach

Problem Summary

Problem 883 asks for the value $T(R_2)$ for remarkable triangles at $R_2=2{,}000{,}000$. The implementations reduce the geometry to counting points in the triangular lattice, where a vector written in basis coordinates $(u,v)$ has squared length

$$Q(u,v)=u^2+uv+v^2.$$

The final answer is not obtained by one direct geometric count: it is a weighted sum over many scaled lattice counts, with Möbius inversion removing non-primitive overcounting.

Mathematical Approach

The triangular-lattice geometry behind remarkable triangles is encoded by the quadratic form $Q(u,v)=u^2+uv+v^2$. The problem therefore becomes an arithmetic-geometric counting problem built from this norm form.

Step 1: Convert the Geometry into Lattice Counts

For a threshold $M$, define

$$A(M)=\#\{(u,v)\in \mathbb{Z}^2:Q(u,v)\le M\},$$

$$A_{\equiv}(M)=\#\{(u,v)\in \mathbb{Z}^2:Q(u,v)\le M,\ u\equiv v\pmod{3}\}.$$

The unrestricted count $A(M)$ measures all admissible lattice vectors up to squared length $M$. The restricted count $A_{\equiv}(M)$ keeps only the congruence class needed when a factor of $3$ has to be absorbed separately. The implementation evaluates these two families for many values of

$$M=\left\lfloor\frac{R_2^2}{t^2}\right\rfloor.$$

Step 2: Count Integer Pairs Efficiently for One Threshold

The key identity is

$$4Q(u,v)=(2v+u)^2+3u^2.$$

After fixing $u$, the condition $Q(u,v)\le M$ becomes

$$ (2v+u)^2\le 4M-3u^2. $$

So $u$ only needs to range over

$$0\le u\le \left\lfloor\sqrt{\frac{4M}{3}}\right\rfloor.$$

If we set

$$s=\left\lfloor\sqrt{4M-3u^2}\right\rfloor,$$

then the admissible integers $v$ satisfy

$$\left\lceil\frac{-u-s}{2}\right\rceil \le v \le \left\lfloor\frac{-u+s}{2}\right\rfloor.$$

The interval length gives the number of points for this $u$. Because replacing $(u,v)$ by $(-u,-v)$ preserves $Q$, the implementation loops over $u\ge0$ and doubles the contribution whenever $u>0$.

Step 3: Explain the Special Congruence Channel

The restricted counter comes from the identity

$$Q(u,v)=u^2+uv+v^2=(u-v)^2+3uv,$$

which implies

$$Q(u,v)\equiv (u-v)^2 \pmod{3}.$$

Therefore

$$3\mid Q(u,v)\iff u\equiv v\pmod{3}.$$

This is why the second channel counts only those values of $v$ in the interval above that satisfy $v\equiv u\pmod{3}$. In Eisenstein-lattice language, the prime $3$ is ramified, so multiples of $3$ must be handled by a separate congruence condition rather than by the ordinary channel.

Step 4: Build the Arithmetic Weights from Prime Factorization

Let

$$n=3^{e_3}\prod_{p\equiv1\pmod{3}}p^{\alpha_p}\prod_{q\equiv2\pmod{3}}q^{\beta_q}.$$

The arithmetic part of the solution distinguishes the three Eisenstein prime types: split primes $p\equiv1\pmod{3}$, inert primes $q\equiv2\pmod{3}$, and the ramified prime $3$. Define

$$\tau_2(n)=\prod_{p\equiv1\pmod{3}}(2\alpha_p+1)\prod_{q\equiv2\pmod{3}}(2\beta_q+1),$$

$$P_1(n)=\prod_{p\equiv1\pmod{3}}(2\alpha_p+1),\qquad \chi(n)=(-1)^{\sum_{q\equiv2\pmod{3}}\beta_q}.$$

The raw multiplicative weight used by the implementation is

$$f(n)= \begin{cases} \tau_2(n)\,(2e_3-1), & e_3\ge 1,\\[4pt] \dfrac{\tau_2(n)+\chi(n)P_1(n)}{2}, & e_3=0. \end{cases}$$

This formula comes directly from how the relevant similarity classes split according to the behavior of primes modulo $3$. The factor $3$ contributes differently from the other primes, which is why it is peeled off first.

Step 5: Use Möbius Inversion to Keep Only Primitive Objects

The raw weight $f(n)$ still counts configurations that share a nontrivial common divisor in the Eisenstein sense. To isolate the primitive contribution, the implementation applies Möbius inversion:

$$h(n)=(f*\mu)(n)=\sum_{d\mid n} f(d)\,\mu\!\left(\frac{n}{d}\right).$$

This Dirichlet convolution removes the contributions that come from scaling up a smaller primitive configuration.

Step 6: Combine the Two Sampling Channels

For each $n\le 3R_2$, the implementations attach a lattice count $F_n$ to the weight $h(n)$:

$$F_n= \begin{cases} A\!\left(\left\lfloor\frac{R_2^2}{n^2}\right\rfloor\right)-1, & 3\nmid n,\\[6pt] A_{\equiv}\!\left(\left\lfloor\frac{R_2^2}{(n/3)^2}\right\rfloor\right)-1, & 3\mid n. \end{cases}$$

The subtraction of $1$ removes the origin. The main accumulated sum is

$$S(R_2)=\sum_{n=1}^{3R_2} h(n)\,F_n.$$

A final symmetry correction is then applied:

$$E(R_2)=\frac{A(R_2^2)-1}{3},\qquad \boxed{T(R_2)=S(R_2)-2E(R_2).}$$

The division by $3$ is valid for the unrestricted nonzero lattice count and reflects the threefold rotational symmetry that would otherwise be counted too many times.

Worked Example: $R_2=1$

This is the smallest checkpoint used by the implementations, and it already shows every moving part. First, $Q(u,v)\le1$ gives the origin plus the six unit directions, so

$$A(1)=7,\qquad A_{\equiv}(1)=1.$$

Hence

$$F_1=A(1)-1=6,\qquad F_3=A_{\equiv}(1)-1=0.$$

Next, compute the raw weights:

$$f(1)=1,\qquad f(2)=\frac{3+(-1)\cdot1}{2}=1,\qquad f(3)=1.$$

Using $\mu(1)=1$, $\mu(2)=-1$, and $\mu(3)=-1$, Möbius inversion gives

$$h(1)=1,\qquad h(2)=f(1)\mu(2)+f(2)\mu(1)=0,\qquad h(3)=f(1)\mu(3)+f(3)\mu(1)=0.$$

Therefore only $n=1$ contributes:

$$S(1)=1\cdot 6=6,\qquad E(1)=\frac{7-1}{3}=2,$$

so the final value is

$$T(1)=6-2\cdot2=2,$$

matching the validation output.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They first build a smallest-prime-factor table and the Möbius function up to $3R_2$. That data is enough to factor every integer $n$, evaluate the multiplicative weight $f(n)$, and then form the primitive weight $h(n)$ by Dirichlet convolution.

Separately, the implementation precomputes the geometric tables for every $1\le t\le R_2$: the unrestricted lattice count and the congruence-restricted lattice count at $M=\left\lfloor R_2^2/t^2\right\rfloor$. The C++ version parallelizes this heavy phase; the Java version performs the same logic directly; the Python version delegates to the compiled computation so that the arithmetic and geometric steps stay identical. Once those tables exist, the program performs one final pass over $n\le 3R_2$, chooses the correct channel according to divisibility by $3$, accumulates $h(n)F_n$, and finally applies the symmetry correction $2E(R_2)$.

Complexity Analysis

For one threshold $M$, the lattice counter runs through $u=0,1,\dots,\left\lfloor\sqrt{4M/3}\right\rfloor$, so a single evaluation costs $O(\sqrt{M})$ time and $O(1)$ extra memory. Summed over all $t\le R_2$, this becomes

$$\sum_{t=1}^{R_2} O\!\left(\sqrt{\frac{R_2^2}{t^2}}\right)=\sum_{t=1}^{R_2} O\!\left(\frac{R_2}{t}\right)=O(R_2\log R_2).$$

The sieve stage is linear or near-linear in $3R_2$, and the Dirichlet-convolution pass costs $O(R_2\log R_2)$ overall. Memory usage is $O(R_2)$: the arithmetic arrays have length $3R_2+1$, and the geometric tables have length $R_2+1$. Parallel execution improves wall-clock time for the geometric phase but does not change the asymptotic bounds.

Footnotes and References

Problem page: https://projecteuler.net/problem=883
Eisenstein integers: Wikipedia — Eisenstein integer
Triangular lattice: Wikipedia — Triangular lattice
Möbius inversion formula: Wikipedia — Möbius inversion formula
Dirichlet convolution: Wikipedia — Dirichlet convolution

Problem 883 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <numeric>
#include <stdexcept>
#include <thread>
#include <utility>
#include <vector>

using namespace std;

namespace {

// floor(sqrt(x)) for 0 <= x <= 2^64-1.
uint64_t isqrt_u64(uint64_t x) {
    long double r = sqrtl(static_cast<long double>(x));
    uint64_t y = static_cast<uint64_t>(r);
    while (static_cast<__int128>(y + 1) * static_cast<__int128>(y + 1) <= x) {
        ++y;
    }
    while (static_cast<__int128>(y) * static_cast<__int128>(y) > x) {
        --y;
    }
    return y;
}

int64_t floor_div2(int64_t a) {
    if (a >= 0) return a >> 1;
    return -(((-a) + 1) >> 1);
}

int64_t ceil_div2(int64_t a) {
    if (a >= 0) return (a + 1) >> 1;
    return -(((-a)) >> 1);
}

int64_t mod_pos(int64_t a, int64_t m) {
    int64_t r = a % m;
    if (r < 0) r += m;
    return r;
}

// Count lattice points (u, v) with u^2 + u v + v^2 <= M.
// Returns (A, Aeq), where A counts all points (including origin),
// and Aeq counts points with u ≡ v (mod 3) (including origin).
pair<uint64_t, uint64_t> count_both(uint64_t M) {
    const uint64_t fourM = 4ULL * M;
    const uint64_t umax = isqrt_u64(fourM / 3ULL);

    uint64_t all = 0;
    uint64_t eq = 0;

    for (uint64_t u = 0; u <= umax; ++u) {
        const uint64_t D = fourM - 3ULL * u * u;
        const uint64_t s = isqrt_u64(D);

        const int64_t uu = static_cast<int64_t>(u);
        const int64_t ss = static_cast<int64_t>(s);
        const int64_t vmin = ceil_div2(-uu - ss);
        const int64_t vmax = floor_div2(-uu + ss);

        uint64_t len = 0;
        if (vmax >= vmin) len = static_cast<uint64_t>(vmax - vmin + 1);

        uint64_t cong = 0;
        if (len) {
            const int r = static_cast<int>(u % 3ULL);
            const int64_t first = vmin + mod_pos(static_cast<int64_t>(r) - vmin, 3);
            if (first <= vmax) {
                cong = static_cast<uint64_t>((vmax - first) / 3 + 1);
            }
        }

        if (u == 0) {
            all += len;
            eq += cong;
        } else {
            all += 2ULL * len;
            eq += 2ULL * cong;
        }
    }
    return {all, eq};
}

// Linear sieve for smallest prime factor and Mobius function.
void linear_sieve_spf_mu(int nmax, vector<int>& spf, vector<int8_t>& mu) {
    spf.assign(nmax + 1, 0);
    mu.assign(nmax + 1, 0);
    vector<int> primes;
    primes.reserve(nmax / 10);

    spf[1] = 1;
    mu[1] = 1;

    for (int i = 2; i <= nmax; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            long long v = 1LL * p * i;
            if (v > nmax) break;
            spf[static_cast<int>(v)] = p;
            if (i % p == 0) {
                mu[static_cast<int>(v)] = 0;
                break;
            }
            mu[static_cast<int>(v)] = static_cast<int8_t>(-mu[i]);
        }
    }
}

// Compute f(n) used in the convolution.
int32_t compute_f_single(int n, const vector<int>& spf) {
    int x = n;
    int e3 = 0;
    while (x % 3 == 0) {
        x /= 3;
        ++e3;
    }

    long long tau_h2 = 1;
    long long prod1 = 1;
    int chi = 1;

    while (x > 1) {
        const int p = spf[x];
        int e = 0;
        while (x % p == 0) {
            x /= p;
            ++e;
        }
        const long long factor = 2LL * e + 1;
        tau_h2 *= factor;

        const int pm3 = p % 3;
        if (pm3 == 1) {
            prod1 *= factor;
        } else if (pm3 == 2) {
            if (e & 1) chi = -chi;
        }
    }

    long long res;
    if (e3 > 0) {
        res = tau_h2 * (2LL * e3 - 1);
    } else {
        res = (tau_h2 + 1LL * chi * prod1) / 2LL;
    }
    return static_cast<int32_t>(res);
}

struct Result {
    int64_t T = 0;
    int64_t S = 0;
    uint64_t points = 0;
};

Result solve_R2(int64_t R2, unsigned numThreads, bool doCheckpoints) {
    const uint64_t N = static_cast<uint64_t>(R2) * static_cast<uint64_t>(R2);
    const int nmax = static_cast<int>(3 * R2);

    vector<int> spf;
    vector<int8_t> mu;
    linear_sieve_spf_mu(nmax, spf, mu);

    vector<int32_t> f(nmax + 1, 0);
    for (int n = 1; n <= nmax; ++n) {
        f[n] = compute_f_single(n, spf);
    }

    // h = f * mu (Dirichlet convolution).
    vector<int32_t> h(nmax + 1, 0);
    for (int k = 1; k <= nmax; ++k) {
        const int8_t muk = mu[k];
        if (muk == 0) continue;
        for (int d = 1, m = k; m <= nmax; ++d, m += k) {
            h[m] += static_cast<int32_t>(muk * f[d]);
        }
    }

    vector<uint64_t> F_base(R2 + 1, 0);
    vector<uint64_t> F_eq(R2 + 1, 0);

    atomic<int> nextT(1);
    const int chunk = 64;

    auto worker = [&]() {
        while (true) {
            const int start = nextT.fetch_add(chunk);
            if (start > R2) break;
            const int end = min<int>(R2, start + chunk - 1);
            for (int t = start; t <= end; ++t) {
                const uint64_t denom = static_cast<uint64_t>(t) * static_cast<uint64_t>(t);
                const uint64_t M = N / denom;
                auto [all, eq] = count_both(M);
                F_base[t] = all - 1ULL;
                F_eq[t] = eq - 1ULL;
            }
        }
    };

    if (numThreads == 0) numThreads = 1;
    vector<thread> threads;
    threads.reserve(numThreads);
    for (unsigned i = 0; i < numThreads; ++i) {
        threads.emplace_back(worker);
    }
    for (auto& th : threads) {
        th.join();
    }

    const uint64_t totalPoints = F_base[1];

    if (doCheckpoints) {
        if (totalPoints % 3ULL != 0ULL) {
            throw runtime_error("Validation failed: totalPoints not divisible by 3");
        }
        for (int t : {1, 2, 3, 10, static_cast<int>(R2)}) {
            if (t <= 0 || t > R2) continue;
            if (F_eq[t] > F_base[t]) {
                throw runtime_error("Validation failed: F_eq[t] > F_base[t]");
            }
        }
    }

    const uint64_t E = totalPoints / 3ULL;

    __int128 S = 0;
    for (int n = 1; n <= nmax; ++n) {
        const int32_t hn = h[n];
        if (hn == 0) continue;

        uint64_t Fn = 0;
        if (n % 3 != 0) {
            if (n <= R2) Fn = F_base[n];
        } else {
            const int t = n / 3;
            Fn = F_eq[t];
        }
        S += static_cast<__int128>(hn) * static_cast<__int128>(Fn);
    }

    const __int128 T = S - 2 * static_cast<__int128>(E);

    Result res;
    res.T = static_cast<int64_t>(T);
    res.S = static_cast<int64_t>(S);
    res.points = totalPoints;
    return res;
}

void run_small_validations(unsigned threads) {
    struct Case { int R2; int64_t expected; };
    const vector<Case> cases = {
        {1, 2},
        {4, 44},
        {20, 1302},
    };

    for (const auto& c : cases) {
        const Result r = solve_R2(c.R2, threads, true);
        if (r.T != c.expected) {
            cerr << "Validation failed for R2=" << c.R2
                 << " (r=" << (c.R2 / 2.0) << "): got " << r.T
                 << ", expected " << c.expected << "\n";
            exit(1);
        }
    }
}

} // namespace

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int64_t R2 = 2'000'000;
    unsigned threads = thread::hardware_concurrency();
    bool validate = true;

    // Optional CLI: ./a.out [R2] [threads] [validate(0/1)]
    if (argc >= 2) R2 = stoll(argv[1]);
    if (argc >= 3) threads = static_cast<unsigned>(stoul(argv[2]));
    if (argc >= 4) validate = (stoi(argv[3]) != 0);

    if (validate) {
        run_small_validations(min(threads, 4u));
    }

    const Result ans = solve_R2(R2, threads, true);
    cout << ans.T << "\n";
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

public class Euler883 {

    static long isqrtU64(long x) {
        long y = (long) Math.sqrt(x);
        while (new java.math.BigInteger(String.valueOf(y + 1)).pow(2)
                .compareTo(new java.math.BigInteger(String.valueOf(x))) <= 0) {
            ++y;
        }
        while (new java.math.BigInteger(String.valueOf(y)).pow(2)
                .compareTo(new java.math.BigInteger(String.valueOf(x))) > 0) {
            --y;
        }
        return y;
    }

    static long floorDiv2(long a) {
        if (a >= 0)
            return a >> 1;
        return -(((-a) + 1) >> 1);
    }

    static long ceilDiv2(long a) {
        if (a >= 0)
            return (a + 1) >> 1;
        return -((-a) >> 1);
    }

    static long modPos(long a, long m) {
        long r = a % m;
        if (r < 0)
            r += m;
        return r;
    }

    static class CountResult {
        long all;
        long eq;

        CountResult(long all, long eq) {
            this.all = all;
            this.eq = eq;
        }
    }

    static CountResult countBoth(long M) {
        long fourM = 4L * M;
        long umax = isqrtU64(fourM / 3L);

        long all = 0;
        long eq = 0;

        for (long u = 0; u <= umax; ++u) {
            long D = fourM - 3L * u * u;
            long s = isqrtU64(D);

            long uu = u;
            long ss = s;
            long vmin = ceilDiv2(-uu - ss);
            long vmax = floorDiv2(-uu + ss);

            long len = 0;
            if (vmax >= vmin)
                len = vmax - vmin + 1;

            long cong = 0;
            if (len > 0) {
                int r = (int) (u % 3L);
                long first = vmin + modPos(r - vmin, 3);
                if (first <= vmax) {
                    cong = (vmax - first) / 3 + 1;
                }
            }

            if (u == 0) {
                all += len;
                eq += cong;
            } else {
                all += 2L * len;
                eq += 2L * cong;
            }
        }
        return new CountResult(all, eq);
    }

    static class SieveResult {
        int[] spf;
        byte[] mu;

        SieveResult(int[] spf, byte[] mu) {
            this.spf = spf;
            this.mu = mu;
        }
    }

    static SieveResult linearSieveSpfMu(int nmax) {
        int[] spf = new int[nmax + 1];
        byte[] mu = new byte[nmax + 1];
        int[] primes = new int[nmax / 10 + 1000];
        int numPrimes = 0;

        spf[1] = 1;
        mu[1] = 1;

        for (int i = 2; i <= nmax; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                if (numPrimes == primes.length) {
                    int[] newPrimes = new int[primes.length * 2];
                    System.arraycopy(primes, 0, newPrimes, 0, primes.length);
                    primes = newPrimes;
                }
                primes[numPrimes++] = i;
                mu[i] = -1;
            }
            for (int j = 0; j < numPrimes; ++j) {
                int p = primes[j];
                long v = 1L * p * i;
                if (v > nmax)
                    break;
                spf[(int) v] = p;
                if (i % p == 0) {
                    mu[(int) v] = 0;
                    break;
                }
                mu[(int) v] = (byte) -mu[i];
            }
        }
        return new SieveResult(spf, mu);
    }

    static int computeFSingle(int n, int[] spf) {
        int x = n;
        int e3 = 0;
        while (x % 3 == 0) {
            x /= 3;
            ++e3;
        }

        long tauH2 = 1;
        long prod1 = 1;
        int chi = 1;

        while (x > 1) {
            int p = spf[x];
            int e = 0;
            while (x % p == 0) {
                x /= p;
                ++e;
            }
            long factor = 2L * e + 1;
            tauH2 *= factor;

            int pm3 = p % 3;
            if (pm3 == 1) {
                prod1 *= factor;
            } else if (pm3 == 2) {
                if ((e & 1) == 1)
                    chi = -chi;
            }
        }

        long res;
        if (e3 > 0) {
            res = tauH2 * (2L * e3 - 1);
        } else {
            res = (tauH2 + 1L * chi * prod1) / 2L;
        }
        return (int) res;
    }

    static long solveR2(long R2) {
        long N = R2 * R2;
        int nmax = (int) (3 * R2);

        SieveResult sieve = linearSieveSpfMu(nmax);
        int[] spf = sieve.spf;
        byte[] mu = sieve.mu;

        int[] f = new int[nmax + 1];
        for (int n = 1; n <= nmax; ++n) {
            f[n] = computeFSingle(n, spf);
        }

        int[] h = new int[nmax + 1];
        for (int k = 1; k <= nmax; ++k) {
            byte muk = mu[k];
            if (muk == 0)
                continue;
            for (int d = 1, m = k; m <= nmax; ++d, m += k) {
                h[m] += muk * f[d];
            }
        }

        long[] FBase = new long[(int) (R2 + 1)];
        long[] FEq = new long[(int) (R2 + 1)];

        // Single threaded for simplicity in translation unless speed is strictly
        // required.
        // 2,000,000 shouldn't be too slow in Java.
        for (int t = 1; t <= R2; ++t) {
            long denom = (long) t * t;
            long M = N / denom;
            CountResult cr = countBoth(M);
            FBase[t] = cr.all - 1L;
            FEq[t] = cr.eq - 1L;
        }

        long totalPoints = FBase[1];
        long E = totalPoints / 3L;

        java.math.BigInteger S = java.math.BigInteger.ZERO;
        for (int n = 1; n <= nmax; ++n) {
            int hn = h[n];
            if (hn == 0)
                continue;

            long Fn = 0;
            if (n % 3 != 0) {
                if (n <= R2)
                    Fn = FBase[n];
            } else {
                int t = n / 3;
                Fn = FEq[t];
            }
            S = S.add(java.math.BigInteger.valueOf(hn).multiply(java.math.BigInteger.valueOf(Fn)));
        }

        java.math.BigInteger T = S.subtract(java.math.BigInteger.valueOf(2 * E));
        return T.longValue();
    }

    public static String solve() {
        return Long.toString(solveR2(2000000L));
    }

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