Problem 386: Maximum Length of an Antichain

View on Project Euler

Project Euler Problem 386 Solution

EulerSolve provides an optimized solution for Project Euler Problem 386, Maximum Length of an Antichain, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each positive integer \(n\), consider the poset of divisors of \(n\) ordered by divisibility. Let \(f(n)\) be the maximum possible size of an antichain in that poset. The implementations in this repository compute $$S(N)=\sum_{n=1}^{N} f(n)$$ for \(N=10^8\). The key simplification is that \(f(n)\) can be read from the middle rank of the divisor lattice, and those middle-rank divisors can be counted by grouping integers according to \(\Omega(n)\), the total number of prime factors with multiplicity. Mathematical Approach Step 1: Rank the divisor lattice by \(\Omega\) Write the prime factorization of \(n\) as $$n=\prod_{i=1}^{r} p_i^{e_i},\qquad \Omega(n)=e_1+\cdots+e_r.$$ Every divisor of \(n\) has the form $$d=\prod_{i=1}^{r} p_i^{a_i},\qquad 0\le a_i\le e_i.$$ If we rank divisors by the sum of their exponents, then the rank of \(d\) is exactly $$\operatorname{rank}(d)=a_1+\cdots+a_r=\Omega(d).$$ Therefore the number of divisors of rank \(m\) is the coefficient of \(x^m\) in the generating polynomial $$P_n(x)=\prod_{i=1}^{r}(1+x+\cdots+x^{e_i}).$$ Step 2: The largest antichain is a middle rank Each factor \(1+x+\cdots+x^{e_i}\) is symmetric and unimodal, and the product of symmetric unimodal polynomials is again symmetric and unimodal. So the coefficient sequence of \(P_n(x)\) increases up to the middle and then decreases symmetrically....

Detailed mathematical approach

Problem Summary

For an integer \(n\), let \(S(n)\) be the set of factors of \(n\), that is, its positive divisors. A subset \(A\subseteq S(n)\) is an antichain if \(A\) has only one element, or if no element of \(A\) divides another. Let \(N(n)\) be the maximum length of an antichain of \(S(n)\). We must find \(\sum N(n)\) over \(1\le n\le 10^8\). [1]

The statement gives one example: \(S(30)=\{1,2,3,5,6,10,15,30\}\); here \(\{2,5,6\}\) is not an antichain, because \(2\) divides \(6\), while \(\{2,3,5\}\) is one. No numeric checkpoint is supplied. The one-element clause makes explicit that a single divisor always counts, so \(N(n)\ge 1\) always, and \(N(1)=1\).

The whole solution in one sentence: a largest antichain of \(S(n)\) consists of the divisors with exactly \(\lfloor\Omega(n)/2\rfloor\) prime factors counted with multiplicity, so pairing each such divisor \(d\) with its cofactor \(n/d\) turns the whole sum into a count of ordered pairs \((a,b)\) with \(ab\le 10^8\) in which \(b\) has as many prime factors as \(a\), or one more — and those pairs are counted by a single monotone sweep.

Mathematical Approach

1. The divisors of n form a grid

Write the prime factorisation \(n=p_1^{e_1}\cdots p_r^{e_r}\). Every divisor is \(d=p_1^{a_1}\cdots p_r^{a_r}\) with \(0\le a_i\le e_i\), and different exponent vectors give different divisors by unique factorisation. So \(S(n)\) is in bijection with the box of vectors \((a_1,\dots,a_r)\) in that range, and \(d\) divides \(d'\) exactly when \(a_i\le a_i'\) for every \(i\): divisibility is the coordinatewise order on a grid.

Let \(\Omega(m)\) count the prime factors of \(m\) with multiplicity, so \(\Omega(d)=a_1+\cdots+a_r\). This function is completely additive: \(\Omega(uv)=\Omega(u)+\Omega(v)\) for all \(u,v\). [6] We call \(\Omega(d)\) the rank of \(d\), and write \(t=\Omega(n)\) for the top rank. For \(n=36=2^2\cdot3^2\) the grid is \(3\times3\), its nine divisors at ranks \(0\) to \(4\); for \(n=30=2\cdot3\cdot5\) it is a \(2\times2\times2\) cube.

2. Divisors of equal rank never divide one another

Suppose \(d\) divides \(d'\) and \(\Omega(d)=\Omega(d')\). Then \(d'=dc\), and additivity gives \(\Omega(c)=0\), so \(c=1\) and \(d'=d\). Hence every rank is an antichain: no two distinct divisors of one rank are comparable.

How big is a rank? Choosing each \(a_i\) independently, the number of divisors of rank \(m\) is the coefficient of \(x^m\) in

$$P_n(x)=\prod_{i=1}^{r}\bigl(1+x+\cdots+x^{e_i}\bigr).$$

For \(n=36\) we get \(P_{36}(x)=(1+x+x^2)^2\), that is \(1+2x+3x^2+2x^3+x^4\): the ranks have sizes \(1,2,3,2,1\), and rank \(2\) is \(\{4,6,9\}\), an antichain of three divisors. The map \(d\mapsto n/d\) is a bijection sending rank \(m\) to rank \(t-m\), so the coefficient list of \(P_n\) always reads the same backwards. Section 6 shows that its largest entry sits in the middle.

3. Chains give upper bounds

A chain is a set of pairwise comparable divisors, written \(d^{(1)}\mid d^{(2)}\mid\cdots\). An antichain and a chain share at most one element, since two elements of a chain are comparable. Hence if \(S(n)\) splits into \(c\) disjoint chains, every antichain has at most \(c\) elements, so \(N(n)\le c\). [2]

Dilworth's theorem says this bound is tight for a suitable splitting, in any finite partial order. [3] We do not need it: for divisors we build a splitting by hand whose chain count an explicit antichain matches, proving equality directly.

4. Symmetric saturated chains

Call a chain in \(S(n)\) saturated if its ranks are consecutive, and symmetric if it starts at some rank \(a\) and ends at rank \(t-a\). Such a chain has \(t-2a+1\) elements, its end ranks adding up to \(t\).

The one fact we need is this. Put \(k=\lfloor t/2\rfloor\). Such a chain, from rank \(a\) to rank \(t-a\), contains exactly one divisor of rank \(k\). Indeed \(a\le t-a\) forces \(a\le t/2\), hence \(a\le k\); and \(k\le t-k\) then gives \(a\le t-k\), that is \(k\le t-a\). So \(k\) lies between the end ranks, and the chain holds one divisor at each rank in that interval.

So if \(S(n)\) splits into symmetric saturated chains, their number equals the number of divisors of rank \(k\). By Section 3 that bounds every antichain, and by Section 2 rank \(k\) is an antichain of that size. Both bounds meet.

5. Building the splitting one prime at a time

We construct one by induction on the number of distinct primes. For \(n=1\) the lone divisor \(1\) is already such a chain, at rank \(0\). For \(n=p^e\) the divisors form the single chain \(1\mid p\mid\cdots\mid p^{e}\), which runs from rank \(0\) to rank \(e=t\) and is saturated and symmetric.

Now let \(n=n'q^{e}\) with \(q\) a prime not dividing \(n'\), and suppose the divisors of \(n'\) are already split into symmetric saturated chains, their top rank being \(t'\). Every divisor of \(n\) is uniquely \(dq^{j}\) with \(d\mid n'\) and \(0\le j\le e\), and \(\Omega(dq^{j})=\Omega(d)+j\). So it suffices to split, for each chain

$$A:\ d_a\mid d_{a+1}\mid\cdots\mid d_{t'-a}$$

of \(n'\), in which \(\Omega(d_i)=i\), the set of products \(d_iq^{j}\). Index these by the cell \((i,j)\) with \(a\le i\le t'-a\) and \(0\le j\le e\), of rank \(i+j\). Write \(L=t'-2a+1\) for the length of \(A\) and \(m=e+1\). For \(s=0,1,\dots,\min(L,m)-1\) define the chain

$$\begin{aligned}C_s:\ &d_aq^{s}\mid d_{a+1}q^{s}\mid\cdots\mid d_{t'-a-s}q^{s}\\&\mid d_{t'-a-s}q^{s+1}\mid\cdots\mid d_{t'-a-s}q^{e}.\end{aligned}$$

In words: multiply \(A\) by \(q^{s}\), walk up it to \(s\) steps before the top, then raise the exponent of \(q\) to the maximum. Each step multiplies by one prime, so \(C_s\) is saturated; it begins at rank \(a+s\) and ends at rank \((t'-a-s)+e\), which add up to \(t'+e\), the top rank of \(n\), so it is symmetric.

They are disjoint and cover the grid. Let \(R_s\) be the cells with \(i\le t'-a-s\) and \(j\ge s\). Then \(R_0\) is the whole grid, \(R_s\) is empty once \(s\ge\min(L,m)\), and removing \(R_{s+1}\) from \(R_s\) leaves exactly the cells of \(R_s\) with \(j=s\) together with those with \(i=t'-a-s\) — precisely \(C_s\). Being the successive differences of a shrinking family, the chains partition \(R_0\).

For an example take \(n'=4\), chain \(1\mid2\mid4\), and adjoin \(q=3\), \(e=2\), reaching \(n=36\). Here \(L=m=3\), giving three chains: \(1\mid2\mid4\mid12\mid36\) at ranks \(0,1,2,3,4\); \(3\mid6\mid18\) at ranks \(1,2,3\); and \(9\) alone at rank \(2\). The nine divisors are used once each, every pair of end ranks adds to \(4\), and each chain holds one divisor of rank \(2\), namely \(4\), \(6\) and \(9\). Those three chains and the antichain \(\{4,6,9\}\) pin \(N(36)=3\) from both sides.

6. The formula for N(n)

Combining Sections 4 and 5: the divisors of any \(n\) split into symmetric saturated chains, so, with \(t=\Omega(n)\),

$$N(n)=\#\{d\mid n:\Omega(d)=\lfloor t/2\rfloor\}.$$

This is the theorem of de Bruijn, van Ebbenhorst Tengbergen and Kruyswijk on the set of divisors of a number. [4] In order-theoretic language, the divisor lattice has the Sperner property: one of its largest antichains consists of elements of a single rank. [5] The same splitting shows the middle rank to be the largest rank, since every rank is an antichain and so has at most as many elements as there are chains.

For odd \(t\) the two middle ranks have equal size by the symmetry of Section 2, so either may be used; the formula fixes the lower one. Small values: \(N(12)=2\) from \(\{2,3\}\), \(N(30)=3\) from \(\{2,3,5\}\) — exactly the antichain the statement exhibits — and \(N(720)=6\). The first twenty values are \(1,1,1,1,1,2,1,1,1,2\) for \(n\le10\), then \(1,2,1,2,2,1,1,2,1,2\), agreeing with the sequence recorded in the OEIS. [7]

7. Each middle divisor is one ordered factor pair

Searching each \(n\)'s divisors separately would be slow. Instead pair every divisor \(d\mid n\) with its cofactor \(b=n/d\): the map \(d\mapsto(d,b)\) is a bijection from \(S(n)\) onto the ordered pairs with \(db=n\), and additivity gives \(\Omega(d)+\Omega(b)=\Omega(n)\).

If \(\Omega(n)=2k\) then \(\lfloor\Omega(n)/2\rfloor=k\), and \(\Omega(d)=k\) says exactly that \(\Omega(d)=\Omega(b)\). If \(\Omega(n)=2k+1\) the floor is again \(k\), and \(\Omega(d)=k\) says exactly that \(\Omega(b)=\Omega(d)+1\). Each \(n\) falls into one case only, so

$$\begin{aligned}N(n)=\ &\#\{(a,b):ab=n,\\ &\quad\Omega(a)=\Omega(b)\}\\ +\ &\#\{(a,b):ab=n,\\ &\quad\Omega(b)=\Omega(a)+1\},\end{aligned}$$

one of the two counts always being zero.

8. One sweep over all pairs

Sum the last identity over \(n\le M\). As a pair \((a,b)\) belongs to exactly one product \(n=ab\), the per-\(n\) searches collapse into two global counts:

$$\begin{aligned}A_M=\ &\#\{(a,b):ab\le M,\\ &\quad\Omega(a)=\Omega(b)\},\end{aligned}$$

$$\begin{aligned}B_M=\ &\#\{(a,b):ab\le M,\\ &\quad\Omega(b)=\Omega(a)+1\},\end{aligned}$$

$$\sum_{n\le M}N(n)=A_M+B_M.$$

Group the integers by their rank: let \(G_k=\{m\le M:\Omega(m)=k\}\), the \(k\)-almost primes up to \(M\). [8] Writing

$$\begin{aligned}H(X,Y)=\ &\#\{(a,b)\in X\times Y:\\ &\quad ab\le M\}\end{aligned}$$

for sets of integers, the two counts become

$$A_M=\sum_{k\ge0}H(G_k,G_k),$$

$$B_M=\sum_{k\ge0}H(G_k,G_{k+1}).$$

Both reduce to one task: count the pairs from two sorted lists whose product is at most \(M\).

9. Worked example: M = 10

Here \(G_0=\{1\}\), \(G_1=\{2,3,5,7\}\), \(G_2=\{4,6,9,10\}\) and \(G_3=\{8\}\). For \(A_{10}\): the pair \((1,1)\); from \(G_1\) the six pairs \((2,2),(2,3),(2,5),(3,2),(3,3),(5,2)\); from \(G_2\) none, since \(4\cdot4=16\) exceeds \(10\); from \(G_3\) none. So \(A_{10}=7\). For \(B_{10}\): the four pairs \((1,2),(1,3),(1,5),(1,7)\), then from \(G_1\times G_2\) only \((2,4)\), because \(2\cdot6=12\) is too large, and nothing from \(G_2\times G_3\). So \(B_{10}=5\), and the sum is \(12\).

Sorting those twelve pairs by product gives one pair each for \(n=1,2,3,4,5,7,8,9\) and two each for \(n=6\) (from \((2,3)\) and \((3,2)\)) and \(n=10\) (from \((2,5)\) and \((5,2)\)) — exactly the values \(N(1),\dots,N(10)\) of Section 6, whose sum is \(12\).

10. Computing the ranks with a sieve

The ranks of all integers up to \(M\) come from one sieve. For \(m\ge2\) let \(\operatorname{spf}(m)\) be the smallest prime factor of \(m\), so \(\operatorname{spf}(p)=p\) for a prime; in the sieve of Eratosthenes the first prime to mark a composite is its smallest prime factor. The loop for \(p\) may start at \(p^2\), since a smaller composite multiple is \(pc\) with \(2\le c\lt p\) and a prime factor of \(c\) marked it already. [9] Additivity then gives the recurrence

$$\begin{aligned}\Omega(1)&=0,\\ \Omega(m)&=\Omega\bigl(m/\operatorname{spf}(m)\bigr)+1,\end{aligned}$$

evaluated in increasing order of \(m\), so the value it refers to is always ready. Appending each \(m\) to its group as \(m\) grows leaves every \(G_k\) sorted for free.

How many groups are there? An integer of rank \(k\) is at least \(2^{k}\). Since \(2^{26}=67108864\) lies below \(10^8\) while \(2^{27}=134217728\) lies above, the largest rank is \(26\); and any rank-\(26\) integer other than \(2^{26}\) would trade a factor \(2\) for a prime at least \(3\), reaching at least \(3\cdot2^{25}=100663296\), already too large. So \(G_{26}=\{2^{26}\}\). The groups are far from uniform: \(|G_1|=\pi(10^8)=5761455\), the largest is \(|G_3|=23727305\), and the sizes then fall away to \(|G_{26}|=1\).

11. Counting pairs with a monotone sweep

Let \(X\) and \(Y\) be sorted lists and count the pairs \((x,y)\in X\times Y\) with \(xy\le M\). For fixed \(x\) the admissible \(y\) are those with \(y\le M/x\), a prefix of \(Y\). If \(x\le x'\) then \(M/x'\le M/x\), so the prefix for \(x'\) lies inside the prefix for \(x\): as \(x\) runs through \(X\) upwards, the index of the last admissible \(y\) only ever moves backwards.

So one pass suffices. Index \(Y\) from zero and start a pointer at \(|Y|-1\); for each \(x\) in increasing order, retreat it while it is nonnegative and the product with that entry exceeds \(M\), then add one more than its value, \(-1\) contributing nothing. With \(M=10\) and \(X=Y=G_1\), the entries \(x=2,3,5,7\) leave the pointer at \(2,1,0,-1\) and contribute \(3,2,1,0\): the six same-rank pairs of Section 9. The pointer retreats at most \(|Y|\) steps in all, so each sweep costs \(O(|X|+|Y|)\) time and \(O(1)\) space, and sweeping \(G_k\) against itself and against \(G_{k+1}\) for every \(k\) produces \(A_M\) and \(B_M\).

How the Code Works

The C++ program [10] reads a struct Options whose limit defaults to \(10^8\); parse_arguments accepts --limit= followed by digits, and --skip-checkpoints. The function gen_primes implements the sieve of Section 10 and returns min_divisor, in which entry \(m\) holds \(\operatorname{spf}(m)\) and min_divisor[1] is \(1\); its name speaks of primes, but the result is the smallest-factor table, and its second parameter is never used.

Inside solve, the byte array degs is filled by the recurrence of Section 10, each integer is appended to by_degs[d], and \(1\) goes into by_degs[0] by hand. The container by_degs starts with \(30\) slots, indexed \(0\) to \(29\), and grows if a rank reaches \(30\), which cannot happen below \(2^{30}\). The counting loop then walks the ranks \(i\) and performs the two sweeps of Section 11 — cur against itself, then cur against nxt — adding high + 1 for every element of cur. Each comparison multiplies by 1LL first, so the product is formed in 64 bits instead of overflowing a 32-bit int; the accumulator ans has type i64.

It is the comparison, not the sum, that truly needs 64 bits: a product is bounded only by \(M^2=10^{16}\), and for the top group C++ really does evaluate \(2^{26}\cdot2^{26}\). The sum is far smaller. Writing \(d(n)\) for the number of divisors, \(N(n)\le d(n)\) and \(\sum_{n\le10^8}d(n)=1857511568\), below \(2^{31}-1\), so a signed 32-bit accumulator would in fact suffice here; i64 is the conservative choice.

The Python version mirrors the same three phases with the lists spf, degs and by_degs, and needs no cast because its integers are unbounded. Its counting loop runs over range(len(by_degs) - 1), so the top occupied group \(G_K\) is never the left-hand list. That is harmless for every \(M\ge2\): every integer of rank \(K\) is at least \(2^{K}\) and every power of two up to \(M\) occurs, so \(2^{K}\le M\lt2^{K+1}\le4^{K}\) once \(K\ge1\), and \(4^{K}\) is the smallest product inside \(G_K\), while \(G_{K+1}\) is empty. At \(M=1\) the loop body never runs and the function returns \(0\) instead of \(1\). Python and Java hard-code their limit, so that observation and the smaller-limit figures below come from separate copies with the constant replaced; C++ returns \(1\) for --limit=1.

The Java version replaces growth by exact allocation. It sieves minDivisor only up to \(\lfloor\sqrt n\rfloor\), then sets minDivisor[i] = i for every entry still zero, marking the primes. It counts the ranks into counts, allocates each byDegs[d] at exactly that length, and fills the arrays in a second pass guided by ptrs. Its degs array is byte[], enough because the largest rank here is \(26\); the sweeps cast one factor to long before multiplying and accumulate into a long.

Only the C++ version carries a self-check: run_checkpoints tests that solve(20) is positive. That confirms the pipeline produces something, not the right thing; solve(20) == 27 would cost nothing. The statement offers no checkpoint, so the small-limit values below were produced independently.

Complexity and Verification

The sieve costs \(O(M\log\log M)\) time. [9] The recurrence for \(\Omega\) and the distribution into groups are one pass each, \(O(M)\). The sweeps make at most \(2M\) outer steps: at \(M=10^8\), C++ makes \(200000000\), Java \(199999999\) and Python \(199999998\), because C++ keeps empty slots past rank \(26\) and visits the top group in both sweeps, Java in the same-rank sweep only, Python in neither. The pointers retreat at most \(|G_k|+|G_{k+1}|\) times per rank, so \(O(M)\) times overall. The whole method runs in \(O(M\log\log M)\) time and \(O(M)\) memory.

Memory binds. In the native arrays of C++ and Java the smallest-factor table needs four bytes per integer, \(381.5\) MiB; the rank array one byte per integer, \(95.4\) MiB; and the groups hold all \(10^8\) integers again as four-byte values, another \(381.5\) MiB. Measured here, C++ finishes in \(0.76\) s at a resident peak of \(1.19\) GiB, the excess being the slack of repeated doubling, which carries \(381.5\) MiB of group data in \(607.5\) MiB of capacity. Java, whose exact arrays avoid that slack, takes \(0.92\) s and \(0.89\) GiB inside a \(900\) MiB heap. Python holds all three structures as lists of references to separate integer objects, so those counts do not describe it: \(31.6\) s and \(5.43\) GiB. A segmented sieve would shrink the sieving workspace, but the groups still hold \(M\) values, so storage stays \(O(M)\). [9]

The following were verified by computation for this article.

  • The formula of Section 6 was not assumed: for every \(n\le1200\) the true width of \(S(n)\) was computed from the definition, as a minimum chain cover found by bipartite matching on the comparability relation, and it equalled the middle-rank count every time.
  • The construction of Section 5 was run for every \(n\) from \(2\) to \(4000\) and checked to partition \(S(n)\) into saturated chains whose end ranks add to \(\Omega(n)\), each holding one divisor of rank \(\lfloor\Omega(n)/2\rfloor\), with as many chains as that rank has elements and that rank the largest.
  • The values \(N(1),\dots,N(20)\) match the OEIS sequence, [7] and summing \(N(n)\) directly, by counting each \(n\)'s divisors of middle rank, gives the same totals as \(A_M+B_M\) for \(M=10,20,100,1000,10^4\); the sweep of Section 11 also agrees with a brute-force double loop over all pairs at \(M=400\).
  • All three programs agree with one another and with those totals — C++ through --limit=, the adapted Python and Java copies through their replaced constant: \(12\), \(27\), \(178\), \(2385\), \(29860\), \(357503\) and \(4153927\) for the limits \(10\), \(20\), \(100\), \(1000\), \(10^4\), \(10^5\) and \(10^6\), and \(47242416\) at \(10^7\) for Java.

Every count is an exact integer, and each integer from \(1\) to \(M\) enters exactly one group, so the group sizes sum to \(M\). The bounds above show that the products and the sum fit the chosen integer types at the required limit; the Python behaviour at a substituted limit of \(1\) is separate from that computation.

Footnotes and References

The references below supply the background facts used above; the chain construction, the pair reduction and the operation counts are derived here.

  1. Project Euler 386 — Maximum Length of an Antichain. The official statement defines S(n) as the set of factors of n, calls a subset an antichain when it has one element or no element divides another, writes N(n) for the maximum length of an antichain, illustrates this with S(30) and the sets {2, 5, 6} and {2, 3, 5}, and asks for the sum of N(n) over 1 ≤ n ≤ 10⁸. No numeric checkpoint is given.
  2. Wikipedia — Antichain. An antichain is a subset of a partial order whose distinct elements are pairwise incomparable, a chain is a totally ordered subset, and the width is the size of the largest antichain. Section 3 needs only that an antichain meets a chain at most once.
  3. Wikipedia — Dilworth's theorem. In any finite partial order the largest antichain has the same size as the smallest partition into chains, the common value being the width. It is context only; Section 5 supplies the equality needed here.
  4. N. G. de Bruijn, C. van Ebbenhorst Tengbergen and D. Kruyswijk, "On the set of divisors of a number", Nieuw Archief voor Wiskunde, series 2, volume 23 (1951), pages 191–193. The original source for Sections 5 and 6: the divisors of an integer split into symmetric chains, and those of middle rank form a largest antichain.
  5. Wikipedia — Sperner's theorem. A graded partial order has the Sperner property when one of its largest antichains consists of elements of a single rank; Sperner's theorem proves this for the subsets of a finite set under inclusion, and Section 6 for divisibility.
  6. Wikipedia — Prime omega function. The function Ω(n) counts the prime factors of n with multiplicity and is completely additive: Ω(uv) = Ω(u) + Ω(v), the identity used in Sections 1, 2, 7 and 10.
  7. OEIS A096825 — maximal size of an antichain in the divisor lattice. The sequence begins 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1, 2, matching Section 6; its comments describe the terms as the divisors whose prime-factor count, with multiplicity, is half that of n, rounded either way.
  8. Wikipedia — Almost prime. An integer is k-almost prime exactly when Ω(n) = k, equivalently when the exponents in its factorisation sum to k; the count up to x is written πₖ(x). These are the groups Gₖ of Section 8.
  9. Wikipedia — Sieve of Eratosthenes. The sieve runs in O(n log log n) operations, the bound following from the reciprocals of the primes; starting the multiples of p at p² is standard; the segmented variant sieves in blocks to reduce memory.
  10. C++ — Euler386.cpp. The linked immutable revision contains gen_primes, solve with the group lists and the two sweeps, and run_checkpoints; it is the implementation record for this article.

Problem 386 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;

struct Options {
    int limit = 100000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(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_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 1;
}

std::vector<int> gen_primes(int n, bool /*minimal_divisors*/) {
    std::vector<int> min_divisor(static_cast<std::size_t>(n + 1), 0);
    min_divisor[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (min_divisor[static_cast<std::size_t>(i)] != 0) {
            continue;
        }
        min_divisor[static_cast<std::size_t>(i)] = i;
        if (1LL * i * i > n) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            if (min_divisor[static_cast<std::size_t>(j)] == 0) {
                min_divisor[static_cast<std::size_t>(j)] = i;
            }
        }
    }
    return min_divisor;
}

i64 solve(int n) {
    const std::vector<int> divs = gen_primes(n, true);
    std::vector<uint8_t> degs(static_cast<std::size_t>(n + 1), 0);
    std::vector<std::vector<int>> by_degs(30);

    by_degs[0].push_back(1);
    for (int i = 2; i <= n; ++i) {
        const int reduced = i / divs[static_cast<std::size_t>(i)];
        const int d = static_cast<int>(degs[static_cast<std::size_t>(reduced)]) + 1;
        degs[static_cast<std::size_t>(i)] = static_cast<uint8_t>(d);
        if (static_cast<int>(by_degs.size()) <= d) {
            by_degs.resize(d + 1);
        }
        by_degs[static_cast<std::size_t>(d)].push_back(i);
    }

    i64 ans = 0;
    for (std::size_t i = 0; i + 1 < by_degs.size(); ++i) {
        const auto& cur = by_degs[static_cast<std::size_t>(i)];

        int low = 0;
        int high = static_cast<int>(cur.size()) - 1;
        while (low < static_cast<int>(cur.size())) {
            while (high >= 0 && 1LL * cur[static_cast<std::size_t>(high)] * cur[static_cast<std::size_t>(low)] > n) {
                --high;
            }
            ans += high + 1LL;
            ++low;
        }

        const auto& nxt = by_degs[static_cast<std::size_t>(i + 1)];
        low = 0;
        high = static_cast<int>(nxt.size()) - 1;
        while (low < static_cast<int>(cur.size())) {
            while (high >= 0 && 1LL * nxt[static_cast<std::size_t>(high)] * cur[static_cast<std::size_t>(low)] > n) {
                --high;
            }
            ans += high + 1LL;
            ++low;
        }
    }

    return ans;
}

bool run_checkpoints() {
    return solve(20) > 0;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        std::cerr << "Checkpoint failed\n";
        return 2;
    }
    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

def solve():
    N = 100_000_000

    # Build min divisor (SPF) sieve
    spf = [0] * (N + 1)
    spf[1] = 1
    for i in range(2, N + 1):
        if spf[i] == 0:
            spf[i] = i
            if i <= N // i:
                for j in range(i*i, N+1, i):
                    if spf[j] == 0:
                        spf[j] = i

    # Compute Omega(n) for each n and group by degree
    degs = [0] * (N + 1)
    by_degs = [[1]]  # degree 0 has {1}
    for i in range(2, N + 1):
        reduced = i // spf[i]
        d = degs[reduced] + 1
        degs[i] = d
        while len(by_degs) <= d:
            by_degs.append([])
        by_degs[d].append(i)

    ans = 0
    for i in range(len(by_degs) - 1):
        cur = by_degs[i]

        # Count pairs within same degree
        low = 0
        high = len(cur) - 1
        while low < len(cur):
            while high >= 0 and cur[high] * cur[low] > N:
                high -= 1
            ans += high + 1
            low += 1

        # Count pairs between consecutive degrees
        nxt = by_degs[i + 1]
        low = 0
        high = len(nxt) - 1
        while low < len(cur):
            while high >= 0 and nxt[high] * cur[low] > N:
                high -= 1
            ans += high + 1
            low += 1

    return str(ans)

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

Java

public class Euler386 {
    static String solve() {
        int n = 100000000;
        int[] minDivisor = new int[n + 1];
        minDivisor[1] = 1;

        int limit = (int) Math.sqrt(n);
        for (int i = 2; i <= limit; i++) {
            if (minDivisor[i] == 0) {
                for (int j = i * i; j <= n; j += i) {
                    if (minDivisor[j] == 0) {
                        minDivisor[j] = i;
                    }
                }
            }
        }
        for (int i = 2; i <= n; i++) {
            if (minDivisor[i] == 0) {
                minDivisor[i] = i;
            }
        }

        byte[] degs = new byte[n + 1];
        int[] counts = new int[32];
        degs[1] = 0;
        counts[0] = 1;

        for (int i = 2; i <= n; i++) {
            int reduced = i / minDivisor[i];
            byte d = (byte) (degs[reduced] + 1);
            degs[i] = d;
            counts[d]++;
        }

        int[][] byDegs = new int[32][];
        for (int i = 0; i < 32; i++) {
            byDegs[i] = new int[counts[i]];
        }

        int[] ptrs = new int[32];
        for (int i = 1; i <= n; i++) {
            int d = degs[i];
            byDegs[d][ptrs[d]++] = i;
        }

        long ans = 0;
        for (int i = 0; i < 31; i++) {
            int[] cur = byDegs[i];
            if (cur == null || cur.length == 0)
                continue;

            int low = 0;
            int high = cur.length - 1;
            while (low < cur.length) {
                while (high >= 0 && (long) cur[high] * cur[low] > n) {
                    high--;
                }
                ans += high + 1;
                low++;
            }

            int[] nxt = byDegs[i + 1];
            if (nxt == null || nxt.length == 0)
                continue;

            low = 0;
            high = nxt.length - 1;
            while (low < cur.length) {
                while (high >= 0 && (long) nxt[high] * cur[low] > n) {
                    high--;
                }
                ans += high + 1;
                low++;
            }
        }

        return Long.toString(ans);
    }

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