Problem 351: Hexagonal Orchards

View on Project Euler

Project Euler Problem 351 Solution

EulerSolve provides an optimized solution for Project Euler Problem 351, Hexagonal Orchards, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary In a hexagonal orchard of order \(n\), every tree lies in one of six congruent triangular sectors around the center. A tree is hidden if another tree lies on the same ray from the center and is closer to the origin. The solution does not compare every pair of trees. Instead it turns visibility into a coprimality condition and then reduces the whole problem to evaluating the summatory Euler totient function. Mathematical Approach Step 1: Reduce the Hexagon to One Sector Because the orchard has sixfold rotational symmetry, it is enough to count hidden trees in one \(60^\circ\) sector and multiply by 6 at the end. In one sector, row \(r\) contains exactly \(r\) trees for \(r=1,2,\dots,n\). Therefore the total number of trees in one sector is the triangular number $$T(n)=\sum_{r=1}^{n} r=\frac{n(n+1)}{2}.$$ The center itself is only the observation point, so it is not included in this count. Step 2: Visibility Becomes a GCD Test Index the trees on row \(r\) by positions \(x=1,2,\dots,r\). In the lattice underlying one hexagonal sector, the tree \((x,r)\) lies on the same ray as \((x/d,r/d)\) whenever \(d\) divides both \(x\) and \(r\). Hence $$\text{tree }(x,r)\text{ is visible from the center} \iff \gcd(x,r)=1.$$ If \(\gcd(x,r)=d>1\), then \((x/d,r/d)\) is a nearer tree on the same line of sight, so \((x,r)\) is hidden....

Detailed mathematical approach

Problem Summary

A hexagonal orchard of order \(n\) is the set of triangular-lattice points inside a regular hexagon of side \(n\), centred at a lattice point and with its six corners at lattice points. An observer stands at the centre, and a non-central point is hidden when another orchard point lies strictly between it and the centre. The statement highlights the thirty hidden points of the orchard of order \(5\), writes \(H(n)\) for their number, gives \(H(5)=30\), \(H(10)=138\), \(H(1\,000)=1177848\), and asks for \(H(10^8)\). [1]

The whole solution in one sentence: in lattice coordinates a point is hidden exactly when its two coordinates have a common divisor greater than \(1\), so of the \(6s\) points on the \(s\)-th hexagonal shell exactly \(6\varphi(s)\) are visible, and the count becomes \(3n(n+1)-6\Phi(n)\) with \(\Phi(n)=\varphi(1)+\cdots+\varphi(n)\) — one number that a recursion over the \(O(\sqrt n)\) distinct quotients \(\lfloor n/k\rfloor\) delivers without listing the \(10^8\) individual totients.

Mathematical Approach

1. Coordinates for the orchard

Put the centre at the origin and take two unit vectors \(60^\circ\) apart,

$$u=(1,0),\qquad v=\left(\tfrac12,\tfrac{\sqrt3}2\right).$$

The triangular lattice is the set of points \(au+bv\) with \(a,b\in\mathbb Z\), written below as the pair \((a,b)\). Both \(u\) and \(v\) point at a corner of the hexagon, and so does \(v-u\), so its corners are the lattice points \((\pm n,0)\), \((0,\pm n)\) and \(\pm(n,-n)\). Its six edges therefore lie on the lines \(a=\pm n\), \(b=\pm n\), \(a+b=\pm n\), and a convex polygon is the intersection of the edge half-planes containing it, so membership means

$$\lVert (a,b)\rVert=\max\bigl(|a|,|b|,|a+b|\bigr)\le n.$$

Call \(\lVert (a,b)\rVert\) the hexagonal distance. Section 3 splits the orchard into shells of \(6s\) points, so its size is

$$1+\sum_{s=1}^{n}6s=3n^2+3n+1,$$

the centred hexagonal number: [5] \(7\), \(19\), \(37\) for \(n=1,2,3\), and \(91\) for \(n=5\), one point being the observer.

2. When is a point hidden?

Hexagonal distance scales like a length: \(\lVert k(a,b)\rVert=k\lVert (a,b)\rVert\) for \(k\ge1\), since each term inside the maximum does. That one fact drives the reduction.

Claim. A point \(P=(a,b)\ne(0,0)\) of the orchard is hidden if and only if \(d=\gcd(|a|,|b|)\ge2\), with \(\gcd(k,0)=k\).

If \(d\ge2\), then \(Q=(a/d,b/d)\) is a lattice point on the segment from the origin to \(P\), strictly nearer, and \(\lVert Q\rVert=\lVert P\rVert/d\le n\), so \(Q\) lies in the orchard and hides \(P\). Conversely let \(d=1\) and let a lattice point \(Q=(c,e)\) lie strictly between the origin and \(P\). Then \(Q=tP\) with \(0\lt t\lt1\), that is \(c=ta\), \(e=tb\). Bézout's identity writes the greatest common divisor of two integers as an integer combination of them — for instance \(1=2\cdot2-3\) — so here \(xa+yb=1\) for some integers \(x,y\), [9] and

$$t=t(xa+yb)=xc+ye\in\mathbb Z,$$

contradicting \(0\lt t\lt1\). The shape of the orchard entered only once: it is convex and centred at the origin, so the blocking point is automatically inside it.

3. Sixfold symmetry, shells and arcs

Rotating by \(60^\circ\) sends \(u\) to \(v\) and \(v\) to \(v-u\), so in coordinates it is

$$\rho(a,b)=(-b,\,a+b),$$

a lattice map with \(\rho^3(a,b)=(-a,-b)\), hence of order \(6\). It permutes \(|a|,|b|,|a+b|\) up to sign, so it preserves hexagonal distance and maps the orchard onto itself, and it preserves \(\gcd(|a|,|b|)\) because \(\gcd(|b|,|a+b|)=\gcd(|a|,|b|)\). Visible points go to visible, hidden to hidden.

The shell \(R_s\) of points with \(\lVert(a,b)\rVert=s\) is the boundary of the hexagon of order \(s\), the union of its six closed edges. The edge from \((s,0)\) to \((0,s)\) lies on \(a+b=s\) with \(a,b\ge0\), so it carries the \(s+1\) lattice points \((s-b,b)\), \(0\le b\le s\), and \(\rho\) carries it to the next edge. Consecutive edges meet in one corner, so

$$|R_s|=6(s+1)-6=6s,$$

and dropping the last corner of each edge leaves six disjoint half-open arcs partitioning \(R_s\), the first being

$$A_s=\{(s-b,b):0\le b\le s-1\}.$$

The shells partition the orchard without its centre, so the rotated arcs do too.

4. One arc counted by Euler's totient

On \(A_s\) the coordinates are \(a=s-b\) and \(b\), and subtracting one argument from the other leaves a greatest common divisor unchanged:

$$\gcd(s-b,\,b)=\gcd(s,\,b).$$

By Section 2 the point \((s-b,b)\) is visible exactly when \(\gcd(s,b)=1\). As \(b\) runs through \(0,1,\dots,s-1\) it runs through a complete set of residues modulo \(s\), so the arc has as many visible points as there are residues coprime to \(s\) — Euler's totient \(\varphi(s)\). [2] For \(s=1\) the only point \((1,0)\) has \(\gcd(1,0)=1\) and is visible, matching \(\varphi(1)=1\).

So \(s-\varphi(s)\) of the \(s\) points are hidden. For \(s=6\) the arc is \((6,0)\), \((5,1)\), \((4,2)\), \((3,3)\), \((2,4)\), \((1,5)\) with \(\gcd(6,b)=6,1,2,3,2,1\); the four hidden ones \((6,0)\), \((4,2)\), \((3,3)\), \((2,4)\) are blocked by \((1,0)\), \((2,1)\), \((1,1)\), \((1,2)\), and \(6-\varphi(6)=4\).

5. The master formula and a worked example

Multiplying by the six rotations and summing over the shells,

$$\begin{aligned}H(n)&=6\sum_{s=1}^{n}\bigl(s-\varphi(s)\bigr)\\&=3n(n+1)-6\Phi(n),\end{aligned}$$

with \(\Phi(n)=\sum_{k\le n}\varphi(k)\) the summatory totient function. For \(n=5\): [4]

  • \(s=1\): \(\varphi(1)=1\), \(0\) per arc, \(0\) on the shell;
  • \(s=2\): \(\varphi(2)=1\), \(1\) per arc, \(6\) on the shell;
  • \(s=3\): \(\varphi(3)=2\), \(1\) per arc, \(6\) on the shell;
  • \(s=4\): \(\varphi(4)=2\), \(2\) per arc, \(12\) on the shell;
  • \(s=5\): \(\varphi(5)=4\), \(1\) per arc, \(6\) on the shell.

Hence \(H(5)=0+6+6+12+6=30\), as the statement says; equivalently \(\Phi(5)=10\) and \(90-60=30\). Likewise \(\Phi(10)=32\) gives \(H(10)=330-192=138\), and \(\Phi(1\,000)=304192\) gives \(H(1\,000)=3003000-1825152\)\({}=1177848\).

6. Gauss's identity turned into a recursion

What remains is one number, \(\Phi(10^8)\). A sieve of all totients would need \(10^8\) array entries and produce \(10^8\) values the formula never uses singly. Start instead from a classical identity of Gauss: for \(m\ge1\),

$$\sum_{d\mid m}\varphi(d)=m,$$

true because reducing the \(m\) fractions \(1/m,\dots,m/m\) to lowest terms produces exactly the denominators \(d\mid m\), each \(\varphi(d)\) times. [2] Sum it over \(m\le n\) and count the pairs \((d,k)\) with \(dk\le n\) in two ways: grouping by \(d\) gives \(\sum_d\varphi(d)\lfloor n/d\rfloor\), grouping by \(k\) gives \(\sum_k\Phi(\lfloor n/k\rfloor)\), since for fixed \(k\) the admissible \(d\) are \(1,\dots,\lfloor n/k\rfloor\). Therefore

$$\sum_{k=1}^{n}\Phi(\lfloor n/k\rfloor)=\frac{n(n+1)}{2},$$

and separating \(k=1\), which is \(\Phi(n)\) itself, leaves a recursion with every argument smaller than \(n\): [3]

$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi(\lfloor n/k\rfloor).$$

For \(n=3\) it reads \(\Phi(3)=6-\Phi(1)-\Phi(1)=4\), which is \(1+1+2\).

7. Equal-quotient blocks and memoisation

As written the sum has \(n-1\) terms, no better than a sieve. The saving lies in the shape of \(\lfloor n/k\rfloor\): it is non-increasing and takes only about \(2\sqrt n\) distinct values, since for \(k\le\sqrt n\) there are at most \(\sqrt n\) indices and for \(k\gt\sqrt n\) the quotient is below \(\sqrt n\). [6] At \(n=10^8\) there are exactly \(19999=2\sqrt n-1\).

Consecutive indices sharing a quotient form a block. If it starts at \(k=\ell\) and \(q=\lfloor n/\ell\rfloor\), it ends at \(r=\lfloor n/q\rfloor\), the largest index whose quotient is still at least \(q\); it contributes \((r-\ell+1)\Phi(q)\), and the next block starts at \(r+1\). For \(n=10\) the loop runs four times:

  • \(\ell=2\): \(q=5\), \(r=2\), subtract \(1\cdot\Phi(5)=10\);
  • \(\ell=3\): \(q=3\), \(r=3\), subtract \(1\cdot\Phi(3)=4\);
  • \(\ell=4\): \(q=2\), \(r=5\), subtract \(2\cdot\Phi(2)=4\);
  • \(\ell=6\): \(q=1\), \(r=10\), subtract \(5\cdot\Phi(1)=5\),

so \(\Phi(10)=55-10-4-4-5=32\), agreeing with Section 5. Caching then makes the recursion cheap, by the nested floor identity

$$\bigl\lfloor\lfloor n/a\rfloor/b\bigr\rfloor=\lfloor n/(ab)\rfloor:$$

writing \(n=ah+r\) and \(h=bj+i\) with \(0\le r\lt a\), \(0\le i\lt b\) gives \(n=abj+(ai+r)\) with \(0\le ai+r\lt ab\), so both sides equal \(j\). Every reachable argument is thus again some \(\lfloor n/k\rfloor\), only those \(\approx2\sqrt n\) values occur, and each is computed once.

Two relatives are worth naming. Let \(\mu(d)\) be the Möbius function: zero when a prime square divides \(d\), else \((-1)^r\) for \(d\) a product of \(r\) distinct primes. Möbius inversion of Gauss's identity gives \(\Phi(n)=\sum_{d\le n}\mu(d)T(\lfloor n/d\rfloor)\) with \(T(m)=m(m+1)/2\), [3] which needs a sieve for \(\mu\) and is no gain here but checks the recursion independently. The partial sums \(M(n)=\sum_{k\le n}\mu(k)\), the Mertens function, obey the same kind of recurrence, with starting value \(1\) in place of \(T(n)\) and the blocks unchanged: [7]

$$M(n)=1-\sum_{k=2}^{n}M(\lfloor n/k\rfloor).$$

8. A sanity check from analysis

The summatory totient satisfies \(\Phi(n)=3n^2/\pi^2+O(n\log n)\), [3] and dividing the master formula by \(n^2\) kills the error term:

$$\begin{aligned}\lim_{n\to\infty}\frac{H(n)}{n^2}&=3-\frac{18}{\pi^2}\\&=1.17621869\ldots,\end{aligned}$$

so the hidden fraction tends to \(1-6/\pi^2\approx0.392\), the complement of the limiting probability \(6/\pi^2\) that two random integers are coprime. The ratio \(H(n)/n^2\) is exactly \(1.177848\) at \(n=1\,000\) and \(1.176874\ldots\) at \(n=3\,000\), enough to catch a lost factor of six.

How the Code Works

The C++ main [8] calls parse_arguments, which reads --order= and --skip-checkpoints and fails with exit code \(1\) on a bad flag; unless the self-test is suppressed it then runs run_checkpoints, and finally prints hidden_points(order).

Sections 6 and 7 live in the class TotientSummatory, whose member memo_ is a std::unordered_map seeded with \(\Phi(0)=0\) and \(\Phi(1)=1\). Its method sum(n) returns a cached value if it has one; otherwise it starts from \(n(n+1)/2\) and walks the blocks — quotient = n / left, right = n / quotient, subtract (right - left + 1) * sum(quotient), jump to left = right + 1. Only that accumulator is an unsigned __int128; the returned and cached values are std::uint64_t, and hidden_points computes in that type. Since \(\Phi(m)\) depends on \(m\) alone, one cache could serve a whole run; the program builds a fresh one in every call of hidden_points.

hidden_points(order) is the master formula verbatim: triangular number, minus summatory.sum(order), times \(6\). brute_hidden_points(order) is the quadratic version of Sections 3 and 4, counting for each \(r\) the \(x\le r\) with \(\gcd(x,r)=1\) by std::gcd. run_checkpoints tests the three published values and the agreement of both methods at order \(200\), and main returns exit code \(2\) if any fails.

Python's solve fixes the order, seeds a local memo dictionary and defines the nested totient_sum over it with the same blocks; unbounded integers make the width question vanish, and the deepest chain of calls, \(\lfloor n/2\rfloor,\lfloor n/4\rfloor,\dots\), descends \(26\) levels at \(n=10^8\), inside the default recursion limit of \(1000\). Java's solve seeds a static HashMap<Long, Long> that outlives the call, and sumTotient mirrors sum in signed long arithmetic. A bound shows this suffices: the terms subtracted for an argument \(m\) are nonnegative and sum to \(T(m)-\Phi(m)\), so each product is at most \(T(m)\le T(n)=5000000050000000\), while the result is below \(3n(n+1)=3\cdot10^{16}\), some \(300\) times under \(2^{63}-1\). Both files print a string from solve and run no checkpoint.

Complexity and Verification

Evaluating one argument \(m\) costs about \(2\sqrt m\) block steps, and the arguments are the \(\approx2\sqrt n\) values \(\lfloor n/k\rfloor\). Splitting them at \(\sqrt n\) and comparing each sum with an integral,

$$\begin{aligned}\sum_{k\le\sqrt n}2\sqrt{n/k}&\approx4n^{3/4},\\ \sum_{m\le\sqrt n}2\sqrt m&\approx\tfrac43n^{3/4},\end{aligned}$$

so the method needs \(\approx\frac{16}{3}n^{3/4}\) arithmetic operations and \(O(\sqrt n)\) stored integers — \(O(n^{3/4})\) time given unit-cost arithmetic and expected constant-time cache lookups. Counting the steps at \(n=10^8\) gives \(5\,254\,261\), split as \(3\,945\,712\) and \(1\,308\,549\) over the two ranges against the predicted \(4\cdot10^6\) and \(1.33\cdot10^6\); the cache ends with \(20\,000\) entries and the recursion descends at most \(26\) levels. The quadratic method would need some \(5\cdot10^{15}\) gcd computations, hence its use only up to order \(200\). On one machine the C++ build (g++ -O2) takes \(0.04\) s including its checkpoints, Java \(0.06\) s and Python \(0.49\) s, all three printing the same value.

Beyond the packaged checkpoints, the following were verified for this article.

  • A brute force using no greatest common divisor — for \(n\le12\) it scans every other point for one collinear with the given one, on the same side of the centre and nearer — finds exactly the hidden points of Section 2, with \(d-1\) blockers each.
  • \(\max(|a|,|b|,|a+b|)\le n\) matches a direct test of the hexagon's edge inequalities for \(|a|,|b|\le n+2\), \(n\le12\); the rotations of \(A_s\) tile \(R_s\) without overlap for \(s\le39\); each shell has \(6s\) points, \(6\varphi(s)\) visible, for \(s\le12\).
  • Enumerating the orchard and testing gcd gives the master formula at \(n=50\), \(137\), \(400\), where \(H(n)=3006\), \(22098\), \(189132\); the first twelve values \(0,6,12,24,30,54\), \(60,84,102,138,144,192\) match the published sequence. [4]
  • Gauss's identity holds for \(m\le299\); the recursion of Section 6 matches a totient sieve for \(n\le800\), the Möbius route for \(n\le2\,000\), and the Mertens recurrence its own sieve for \(n\le2\,000\).

Three limits deserve plain statement. Nothing checks overflow, and the \(128\)-bit accumulator inside sum does not make a large --order= safe: the result leaves the unsigned \(64\)-bit range near \(n\approx3.96\cdot10^9\), and order * (order + 1ULL) at \(n=2^{32}\). Decimal parsing is unchecked, so --order=18446744073709551617 wraps to order \(1\) and prints \(0\). And since the Python and Java files run no checkpoints, a change to either passes unnoticed unless its output is compared with the C++ one.

Footnotes and References

The references supply the background facts; the coordinates, the hidden-point criterion, the shell decomposition and the operation counts are derived here.

  1. Project Euler 351 — Hexagonal Orchards. Defines the orchard of order n as a triangular lattice made up of points within a regular hexagon with side n, highlights those hidden from the centre by a closer point, and gives H(5) = 30, H(10) = 138 and H(1000) = 1177848 before asking for order 100000000.
  2. Wikipedia — Euler's totient function. Gives φ(1) = 1, Gauss's identity that the φ(d) over the divisors of n sum to n with the proof by reducing k/n, and the limiting coprimality probability 6/π².
  3. OEIS A002088 — sum of the totient function. Lists Φ(n) = 0, 1, 2, 4, 6, 10, 12, 18, 22, 28, 32 for n = 0…10, the estimate 3n²/π² + O(n log n), the recursion a(n) = n(n+1)/2 − Σ a(⌊n/k⌋) over k ≥ 2 used here, and the Möbius form ½ Σ μ(k)⌊n/k⌋⌊1 + n/k⌋.
  4. OEIS A216453 — hidden points of a hexagonal orchard. The points hidden from the central point by a closer point in an orchard of order n: it cites this problem, gives a(n) = 6(C(n+1,2) − Σ φ(i)), and lists 57 terms from 0, 6, 12, 24, 30, 54 on, with a(50) = 3006.
  5. OEIS A003215 — centred hexagonal numbers. The values 3n(n+1) + 1 = 1, 7, 19, 37, 61, 91, …, the crystal ball sequence of the hexagonal lattice: the lattice points within distance n of a point.
  6. Wikipedia — divisor summatory function. Names the Dirichlet hyperbola method and computes Σ ⌊x/k⌋ in O(√x) operations by splitting the range at √x — the background for Section 7's split into small indices and small quotients.
  7. Wikipedia — Mertens function. Defines M(n) as the partial sums of the Möbius function, records the companion identity Σ M(n/d) = 1, and reports the O(x^(2/3)(log log x)^(1/3)) method of Deléglise and Rivat, better than the O(n^(3/4)) here at the price of a sieve.
  8. C++ — Euler351.cpp. The linked immutable revision holds TotientSummatory, hidden_points, brute_hidden_points and the checkpoints above; the implementation record for this article.
  9. Wikipedia — Bézout's identity. For integers a, b with greatest common divisor d there are integers x, y with ax + by = d, proved through the least positive value of ax + by; Section 2 uses the coprime case d = 1.

Problem 351 source code

C++

#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <unordered_map>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

struct Options {
    u64 order = 100'000'000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    u64 parsed = 0ULL;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--order=", options.order)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.order >= 1ULL;
}

class TotientSummatory {
public:
    TotientSummatory() {
        memo_[0ULL] = 0ULL;
        memo_[1ULL] = 1ULL;
    }

    u64 sum(const u64 n) {
        const auto it = memo_.find(n);
        if (it != memo_.end()) {
            return it->second;
        }
        u128 total = static_cast<u128>(n) * static_cast<u128>(n + 1ULL) / 2ULL;
        for (u64 left = 2ULL; left <= n;) {
            const u64 quotient = n / left;
            const u64 right = n / quotient;
            total -= static_cast<u128>(right - left + 1ULL) * static_cast<u128>(sum(quotient));
            left = right + 1ULL;
        }
        const u64 result = static_cast<u64>(total);
        memo_[n] = result;
        return result;
    }

private:
    std::unordered_map<u64, u64> memo_;
};

u64 hidden_points(const u64 order) {
    TotientSummatory summatory;
    const u64 triangle = order * (order + 1ULL) / 2ULL;
    const u64 visible_per_sector = summatory.sum(order);
    const u64 hidden_per_sector = triangle - visible_per_sector;
    return hidden_per_sector * 6ULL;
}

u64 brute_hidden_points(const int order) {
    u64 hidden_per_sector = 0ULL;
    for (int r = 1; r <= order; ++r) {
        int visible = 0;
        for (int x = 1; x <= r; ++x) {
            if (std::gcd(x, r) == 1) {
                ++visible;
            }
        }
        hidden_per_sector += static_cast<u64>(r - visible);
    }
    return hidden_per_sector * 6ULL;
}

bool run_checkpoints() {
    if (hidden_points(5ULL) != 30ULL) {
        std::cerr << "Checkpoint failed for H(5)" << '\n';
        return false;
    }
    if (hidden_points(10ULL) != 138ULL) {
        std::cerr << "Checkpoint failed for H(10)" << '\n';
        return false;
    }
    if (hidden_points(1'000ULL) != 1'177'848ULL) {
        std::cerr << "Checkpoint failed for H(1000)" << '\n';
        return false;
    }
    if (hidden_points(200ULL) != brute_hidden_points(200)) {
        std::cerr << "Checkpoint failed for brute-force cross-check at order=200" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << hidden_points(options.order) << '\n';
    return 0;
}

Python

def solve():
    order = 100_000_000

    memo = {0: 0, 1: 1}

    def totient_sum(n):
        if n in memo:
            return memo[n]
        total = n * (n + 1) // 2
        left = 2
        while left <= n:
            q = n // left
            right = n // q
            total -= (right - left + 1) * totient_sum(q)
            left = right + 1
        memo[n] = total
        return total

    triangle = order * (order + 1) // 2
    visible = totient_sum(order)
    hidden = triangle - visible
    return str(hidden * 6)

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

Java

import java.util.HashMap;

public class Euler351 {
    static HashMap<Long, Long> memo = new HashMap<>();

    static long sumTotient(long n) {
        Long cached = memo.get(n);
        if (cached != null)
            return cached;
        long total = n * (n + 1) / 2;
        long left = 2;
        while (left <= n) {
            long quotient = n / left;
            long right = n / quotient;
            total -= (right - left + 1) * sumTotient(quotient);
            left = right + 1;
        }
        memo.put(n, total);
        return total;
    }

    public static String solve() {
        memo.put(0L, 0L);
        memo.put(1L, 1L);
        long order = 100000000L;
        long triangle = order * (order + 1) / 2;
        long visible = sumTotient(order);
        long hidden = triangle - visible;
        long ans = hidden * 6;
        return String.valueOf(ans);
    }

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