Problem 420: $2 \times 2$ Positive Integer Matrix

Detailed mathematical approach

Problem Summary

A positive integer matrix is a matrix all of whose entries are positive integers. Some \(2\times2\) positive integer matrices are the square of a positive integer matrix in two different ways; the statement gives the example

$$\begin{aligned}\begin{pmatrix}40&12\\48&40\end{pmatrix}&=\begin{pmatrix}2&3\\12&2\end{pmatrix}^2\\&=\begin{pmatrix}6&1\\4&6\end{pmatrix}^2.\end{aligned}$$

Let \(F(N)\) be the number of \(2\times2\) positive integer matrices whose trace — the sum of the two entries on the main diagonal — is less than \(N\) and which can be written as the square of a positive integer matrix in two different ways. We are told that \(F(50)=7\) and \(F(1000)=1019\), and we must find \(F(10^7)\). ^[1]

The whole solution in one sentence: a matrix with two positive integer square roots is pinned down by the traces of those two roots, and once those traces are chosen the only freedom left is to split one small integer into two factors, so the whole count turns into a sum of divisor counts over a few million cases.

Mathematical Approach

1. The square of a 2×2 matrix

Let \(B=\begin{pmatrix}a&b\\c&d\end{pmatrix}\) and let \(t=a+d\) be its trace. Its determinant is \(\det B=ad-bc\), and \(I\) is the identity matrix, with ones on the main diagonal and zeros elsewhere. The Cayley–Hamilton theorem says \(B^2-tB+(\det B)I=0\), ^[2] and multiplying out confirms

$$B^2=\begin{pmatrix}a^2+bc&bt\\ct&d^2+bc\end{pmatrix}.$$

Two features of this shape drive everything that follows: both off-diagonal entries of \(B^2\) are multiples of the trace \(t\), and the difference of the diagonal entries is \((a-d)(a+d)\). The trace of the square is \(a^2+d^2+2bc\). For \(B=\begin{pmatrix}2&3\\12&2\end{pmatrix}\) we get \(t=4\), \(a^2+bc=4+36=40\), \(bt=12\) and \(ct=48\), which is the example matrix above.

2. A square root is determined by its trace

Fix a target \(A=\begin{pmatrix}A_{11}&A_{12}\\A_{21}&A_{22}\end{pmatrix}\) and write \(S=A_{11}+A_{22}\) for its trace and \(D=A_{11}-A_{22}\). Suppose \(B^2=A\) with \(B\) as in Section 1 and \(t=a+d\). Comparing entries with Section 1 gives four equations:

• \(b=A_{12}/t\) and \(c=A_{21}/t\);
• \(a-d=D/t\), because \(A_{11}-A_{22}=(a-d)t\);
• \(a+d=t\), so \(a\) and \(d\) are determined as well.

Every entry of \(B\) is therefore a function of \(A\) and of the single number \(t\). Two different square roots of the same matrix must consequently have two different traces.

3. Only two traces are possible

Using \(a^2+d^2=\tfrac12\bigl((a+d)^2+(a-d)^2\bigr)\) and \(bc=A_{12}A_{21}/t^2\), the trace of \(A\) becomes

$$S=\frac{t^2+D^2/t^2}{2}+\frac{2A_{12}A_{21}}{t^2}.$$

Multiplying by \(2t^2\) turns this into a quadratic equation in \(t^2\),

$$t^4-2St^2+D^2+4A_{12}A_{21}=0,$$

and completing the square, with \(S^2-D^2=4A_{11}A_{22}\), gives the compact form

$$(t^2-S)^2=4\det A.$$

This single identity carries three consequences. First, \(\det A=(\det B)^2\) is a perfect square. Second, \(t^2=S\pm2\sqrt{\det A}\) leaves at most two values of \(t^2\), hence at most two positive values of \(t\), hence by Section 2 at most two positive integer square roots: "in two different ways" therefore means exactly two ways, and counting the matrices is the same as counting the pairs of roots. Third, if \(t_1\lt t_2\) are the two traces then

$$t_1^2+t_2^2=2S,\qquad t_2^2-t_1^2=4\sqrt{\det A}.$$

The example has \(S=80\) and \(\det A=40\cdot40-12\cdot48=1024\), so \(t^2\) is \(16\) or \(144\) and the two traces are \(t_1=4\) and \(t_2=12\) — and indeed the two roots given in the statement have traces \(4\) and \(12\). A matrix with \(\det A=0\) admits only one trace and is never counted.

4. The ratio of the two traces

Let \(B_1=\begin{pmatrix}a&b\\c&d\end{pmatrix}\) and \(B_2=\begin{pmatrix}e&f\\g&h\end{pmatrix}\) be the two roots, with traces \(t_1=a+d\lt t_2=e+h\). Put \(n=\gcd(t_1,t_2)\), the greatest common divisor of the two traces. Dividing both traces by \(n\) cancels every common factor, so

$$t_1=qn,\qquad t_2=pn,\qquad \gcd(p,q)=1,\ p\gt q\ge1.$$

The first relation of Section 3 now fixes the trace of \(A\) completely:

$$\operatorname{tr}(A)=\frac{t_1^2+t_2^2}{2}=\frac{(p^2+q^2)n^2}{2}.$$

It depends on \(p\), \(q\) and \(n\) only, so the condition \(\operatorname{tr}(A)\lt N\) is simply \((p^2+q^2)n^2\le 2N-1\). For the example, \(t_1=4\) and \(t_2=12\) give \(n=4\), \(q=1\), \(p=3\), and \((9+1)\cdot16/2=80\) is its trace.

5. The diagonal entries

Write \(r_1=a-d\) and \(r_2=e-h\). Both roots produce the same difference of diagonal entries, so by Section 1

$$A_{11}-A_{22}=r_1t_1=r_2t_2,$$

and dividing by \(n\) gives \(qr_1=pr_2\).

Since \(p\) divides \(qr_1\) and \(\gcd(p,q)=1\), \(p\) divides \(r_1\); ^[3] write \(r_1=pm\) with an integer \(m\), and then \(r_2=qm\). The four diagonal entries are now

$$a=\frac{qn+pm}{2},\qquad d=\frac{qn-pm}{2},$$

$$e=\frac{pn+qm}{2},\qquad h=\frac{pn-qm}{2}.$$

The integer \(m\) may be negative or zero, and its sign matters. Replacing a nonzero \(m\) by \(-m\) exchanges \(a\) with \(d\) and \(e\) with \(h\); by Section 1 the square then has its two diagonal entries exchanged, and since \(a\ne d\) that is a different matrix \(A\).

6. The off-diagonal entries

By Section 1 the off-diagonal entries of \(A\) are \(A_{12}=bt_1=ft_2\), that is \(bq=fp\) after dividing by \(n\). As in Section 5, \(p\) divides \(b\), so \(b=pu\) and \(f=qu\); likewise \(c=pv\) and \(g=qv\). Here \(u\) and \(v\) are positive integers, because \(b\) and \(c\) are positive multiples of \(p\). In terms of the new parameters

$$A_{12}=pqnu,\qquad A_{21}=pqnv.$$

7. The one equation that ties the parameters together

So far only the off-diagonal entries and the difference of the diagonal entries have been matched. The remaining condition is that the upper-left entries agree, \(a^2+bc=e^2+fg\), that is \(a^2+p^2uv=e^2+q^2uv\). From Section 5,

$$\begin{aligned}a^2-e^2&=\frac{(qn+pm)^2-(pn+qm)^2}{4}\\&=\frac{(q^2-p^2)(n^2-m^2)}{4},\end{aligned}$$

because the mixed terms \(2pqnm\) cancel. The entry equation reads \(a^2-e^2=(q^2-p^2)uv\), so dividing by \(q^2-p^2\ne0\) leaves one Diophantine condition,

$$uv=\frac{n^2-m^2}{4}=:\ell.$$

The lower-right entries give nothing new: \((qn-pm)^2-(pn-qm)^2\) equals \((q^2-p^2)(n^2-m^2)\) as well. So a quadruple \((p,q,n,m)\) fixes the number \(\ell\), and the only remaining freedom is the choice of an ordered pair \((u,v)\) of positive integers with \(uv=\ell\).

8. Positivity

The diagonal entries of \(B_1\) are positive exactly when \(qn+pm\) and \(qn-pm\) are both positive, that is when

$$p|m|\lt qn,$$

which is the same as \(|m|\le\lfloor(qn-1)/p\rfloor\).

This one inequality also settles the rest. It gives \(p|m|\lt qn\le pn\), hence \(|m|\lt n\), so \(\ell\gt0\) automatically; and \(q|m|\le p|m|\lt qn\le pn\) makes \(e\) and \(h\) positive too. Note that \(n=1\) can never occur: it would force \(m=0\), and for \(m=0\) the parity rule of Section 9 requires \(n\) to be even.

9. Parity

The four diagonal entries are integers exactly when \(qn+pm\) and \(pn+qm\) are even; the other two numerators differ from these by \(2pm\) and \(2qm\). Coprime \(p\) and \(q\) are not both even, which leaves two cases.

• Opposite parity. Say \(q\) is odd and \(p\) even. Then \(qn+pm\) is even only if \(n\) is even, and \(pn+qm\) only if \(m\) is even; the case of odd \(p\) and even \(q\) is symmetric. So \(n\) and \(m\) are both even.
• Both odd. Both conditions reduce to \(n\equiv m\pmod 2\).

In either case \(4\) divides \(n^2-m^2\), so \(\ell\) is a positive integer: for even \(n\) and \(m\) this is clear, and when both are odd, odd squares are congruent to \(1\) modulo \(8\), so \(n^2-m^2\) is even a multiple of \(8\).

10. Why each quadruple contributes \(\tau(\ell)\)

Call a quadruple \((p,q,n,m)\) admissible when \(\gcd(p,q)=1\), \(p\gt q\ge1\), \((p^2+q^2)n^2\le2N-1\), \(p|m|\lt qn\) and the parity rule of Section 9 holds, and let \(\mathcal A\) be the set of admissible quadruples. Sections 4 to 9 show that every matrix counted by \(F(N)\) determines exactly one admissible quadruple together with one ordered factorisation \(uv=\ell\): the two traces give \(n\), \(q\) and \(p\), the smaller root gives \(m=(a-d)/p\), \(u=b/p\) and \(v=c/p\).

Conversely, take any admissible quadruple and any ordered factorisation \(uv=\ell\). The eight entries of Sections 5 and 6 are then positive integers, the traces \(qn\ne pn\) differ, so the two matrices are different, and reading Section 7 backwards shows that their squares agree in every entry. The common square has trace \((p^2+q^2)n^2/2\lt N\) and is therefore counted by \(F(N)\).

The two directions are inverse to each other: a matrix has only the two roots of Section 3, so the recipe recovers the very quadruple and factorisation it was built from. The count is therefore a sum, over admissible quadruples, of the number of ordered factorisations of \(\ell\), which is the number of divisors \(\tau(\ell)=\sigma_0(\ell)\), ^{[4] [5]} since each divisor \(u\) of \(\ell\) gives one pair \((u,\ell/u)\):

$$F(N)=\sum_{(p,q,n,m)\in\mathcal A}\tau\Bigl(\frac{n^2-m^2}{4}\Bigr).$$

11. The bounds needed before the sum

Among coprime pairs with \(p\gt q\ge1\) the smallest value of \(p^2+q^2\) is \(5\), at \((p,q)=(2,1)\). The trace bound therefore limits every parameter:

$$\begin{aligned}n&\le\Bigl\lfloor\sqrt{\tfrac{2N-1}{5}}\Bigr\rfloor,\\ \ell&\le\frac{n^2}{4},\qquad p^2\lt 2N.\end{aligned}$$

For \(N=10^7\) this means \(n\le1999\), \(\ell\le999\,000\) and \(p\le4472\). A table of \(\tau\) up to \(999\,000\), built once, answers every one of the millions of lookups in the sum.

12. Worked example: \(F(50)=7\)

Here \((p^2+q^2)n^2\le99\). Since \(n=1\) never occurs (Section 8) we need \(p^2+q^2\le24\), which leaves the coprime pairs \((2,1)\), \((3,1)\), \((3,2)\) and \((4,1)\) with \(n=2\), plus \((2,1)\) with \(n=4\). In each case the bound \(p|m|\lt qn\) together with the parity rule leaves \(m=0\) alone. The trace bound also allows \((3,1)\) with \(n=3\), but these \(p\) and \(q\) are both odd, so \(m\) would have to be odd with \(3|m|\lt3\), which is impossible. That gives five admissible quadruples:

• \((2,1,2,0)\): \(\ell=1\), \(\tau=1\), trace \(10\);
• \((3,1,2,0)\): \(\ell=1\), \(\tau=1\), trace \(20\);
• \((3,2,2,0)\): \(\ell=1\), \(\tau=1\), trace \(26\);
• \((4,1,2,0)\): \(\ell=1\), \(\tau=1\), trace \(34\);
• \((2,1,4,0)\): \(\ell=4\), \(\tau=3\), trace \(40\).

The total is \(1+1+1+1+3=7=F(50)\). The four quadruples with \(\ell=1\) force \(u=v=1\); for instance \((2,1,2,0)\) gives \(\begin{pmatrix}1&2\\2&1\end{pmatrix}\) and \(\begin{pmatrix}2&1\\1&2\end{pmatrix}\), both squaring to \(\begin{pmatrix}5&4\\4&5\end{pmatrix}\). The quadruple \((2,1,4,0)\) has \(\ell=4\) and \(\tau(4)=3\), so its three ordered factorisations \((u,v)=(1,4),(2,2),(4,1)\) give the three matrices of trace \(40\):

$$\begin{pmatrix}20&8\\32&20\end{pmatrix},\ \begin{pmatrix}20&16\\16&20\end{pmatrix},\ \begin{pmatrix}20&32\\8&20\end{pmatrix},$$

the first of them being \(\begin{pmatrix}2&2\\8&2\end{pmatrix}^2=\begin{pmatrix}4&1\\4&4\end{pmatrix}^2\). The statement's own example is the quadruple \((3,1,4,0)\) with \(\ell=4\) and \((u,v)=(1,4)\); its trace is \(80\), so it is counted by \(F(1000)\) but not by \(F(50)\).

13. A worked example with \(m\ne0\)

Trace \(90\) needs \((p^2+q^2)n^2=180\), whose only coprime solution is \(p=2\), \(q=1\), \(n=6\). These have opposite parity, so \(m\) is even, and \(2|m|\lt6\) leaves \(m\in\{-2,0,2\}\), with \(\ell=(36-m^2)/4\) equal to \(8\), \(9\) and \(8\). Since \(\tau(8)=4\) and \(\tau(9)=3\), exactly \(4+3+4=11\) counted matrices have trace \(90\), which agrees with \(F(91)-F(90)=27-16=11\). The two signs of \(m\) give different matrices, as announced in Section 5: with \(m=-2\) the smaller root is \(\begin{pmatrix}1&2u\\2v&5\end{pmatrix}\) and with \(m=+2\) it is \(\begin{pmatrix}5&2u\\2v&1\end{pmatrix}\). Taking \(m=-2\) and \((u,v)=(1,8)\),

$$\begin{aligned}\begin{pmatrix}33&12\\96&57\end{pmatrix}&=\begin{pmatrix}1&2\\16&5\end{pmatrix}^2\\&=\begin{pmatrix}5&1\\8&7\end{pmatrix}^2.\end{aligned}$$

How the Code Works

All three programs evaluate the sum of Section 10 with the bounds of Section 11, in four nested loops over \(p\), \(q\), \(n\) and \(m\).

In the C++ file, parse_arguments reads --n= (default \(10^7\)) and --skip-checkpoints and insists on \(N\ge2\). isqrt_i64 takes a long double square root and then corrects it upwards and downwards until it is the exact integer square root, so the floors of Section 11 are never off by one. F_value does the work:

• it computes p_max, n_global and l_max from Section 11;
• it fills a std::vector<int> named tau of l_max + 1 entries with the divisor sieve: for every \(d\) it increments tau[x] at \(x=d,2d,3d,\dots\), so afterwards tau[x] is \(\tau(x)\);
• it loops over \(p\) from \(2\) to p_max and \(q\) from \(1\) to \(p-1\), keeping the coprime pairs with std::gcd; ^[6]
• for each pair it sets pq2 \(=p^2+q^2\), n_max \(=\lfloor\sqrt{(2N-1)/pq2}\rfloor\) and the flag mixed_parity \(=\) ((p ^ q) & 1) != 0, which is true exactly when \(p\) and \(q\) have opposite parity;
• for each \(n\) it skips odd \(n\) when mixed_parity holds (Section 9), computes m_max \(=\lfloor(qn-1)/p\rfloor\) (Section 8) and the required parity of \(m\), which is \(n\bmod2\) when \(p\) and \(q\) are both odd and \(0\) otherwise;
• it starts at \(m=-\)m_max, moves one step up if that value has the wrong parity, and then steps by \(2\), computing l_num = n2 - m * m and l = l_num / 4 and adding tau[l] to an unsigned 64-bit accumulator. The loop variable is called nn, with n2 = nn * nn, because the name n is taken by the bound \(N\).

The parity test m & 1 is applied to a negative m, so it relies on two's complement, the only signed representation C++ has allowed since C++20: ^[7] \(-3\) has last bit \(1\) and \(-4\) has last bit \(0\), which is what the correction needs. Two guards in the inner code can never fire, because Section 8 already rules the cases out: m_max is never negative, and \(n^2-m^2\) is never \(\le0\).

brute_count_small is an independent check of the whole derivation for a small bound. It runs over all \(a,d\le\lfloor\sqrt{N-1}\rfloor\), computes how much of the trace budget is left for \(2bc\), enumerates every \(b\) and \(c\) within that budget, builds the square as in Section 1 and packs its four entries into one 64-bit key, sixteen bits each, in an std::unordered_map; the matrices reached at least twice are exactly those with two roots. With \(N=200\) it inserts \(29\,143\) squares, \(29\,065\) of them distinct, finds the \(78\) keys reached twice, and never stores an entry above \(759\). run_checkpoints then verifies \(F(50)=7\) and \(F(1000)=1019\) against the statement and F_value(200) against this brute force, and main prints \(F(N)\) for the requested \(N\). Each checkpoint call rebuilds its own small \(\tau\) table.

The Python file has the same four loops, with math.gcd and math.isqrt; ^[8] its wrapper around math.isqrt with two correction loops is redundant, since math.isqrt is already exact, and its integers are unbounded, so no overflow question arises. The Java file uses a long accumulator, an int[] for tau and its own Euclidean gcd. Its start correction tests (m & 1) != parity && (m & 1) != -parity, which is equivalent to the C++ test, because m & 1 is \(0\) or \(1\) for every Java int. ^[9] Neither port contains the checkpoints or the brute force; both compute \(F(10^7)\) only.

Both the C++ and the Java version keep products such as p2 + q * q in a 32-bit int. At the requested size that is safe, the largest of them being \(39\,988\,625\) at \(p=4472\), \(q=4471\), but parse_arguments sets no upper bound on \(N\): the accepted bound \(536\,903\,681\) would reach the coprime pair \(p=32769\), \(q=32768\), whose \(p^2+q^2=2\,147\,549\,185\) no longer fits such an integer. Python is immune, since its integers grow as needed.

Complexity and Verification

The divisor sieve performs \(\sum_{d\le L}\lfloor L/d\rfloor\) increments, which is \(L\ln L+O(L)\) because the harmonic numbers satisfy \(H_L=\ln L+\gamma+O(1/L)\). ^[10] At \(N=10^7\) the bound is \(L=999\,000\) and the sieve makes \(13\,955\,099\) increments. In C++ and Java the \(L+1\) entries are four-byte integers and occupy \(3.8\) MiB; Python holds tau as a list of references to integer objects, which costs considerably more for the same entries.

The four loops cost, at \(N=10^7\): \(9\,997\,156\) calls of gcd, of which \(4\,774\,679\) pairs are coprime and admit some \(n\); \(7\,850\,162\) iterations of the \(n\) loop, of which \(3\,924\,300\) survive the parity filter; and \(9\,866\,844\) admissible quadruples, one table lookup each.

These counts can be bounded in advance. Fix \(n\). The trace bound forces \(p^2n^2\lt2N\), so at most \(O(N/n^2)\) pairs \((p,q)\) with \(q\lt p\) remain, and for each pair Section 8 leaves at most \(n\) values of \(m\); dropping the coprimality and parity filters only increases both counts. Summing over \(n\le\sqrt{2N}\),

$$\sum_{n\le\sqrt{2N}}\frac{N}{n^2}\cdot n=O(N\log N),$$

and summing \(O(N/n^2)\) alone bounds the \((p,q,n)\) iterations by \(O(N)\), because \(\sum1/n^2\) converges. The two outer loops run over \(\sum_{p\le\sqrt{2N}}(p-1)=O(N)\) pairs, each costing one gcd, which takes \(O(\log N)\) steps. Together with the sieve's \(L\ln L\) and \(L\approx N/10\), the method needs \(O(N\log N)\) operations on numbers of \(O(\log N)\) bits and \(O(N)\) words of storage. The measured quadruple counts \(4949\), \(66\,047\), \(824\,272\) and \(9\,866\,844\) at \(N=10^4,10^5,10^6,10^7\), that is \(0.49N\) to \(0.99N\) with about \(0.16N\) added per decade, behave as that bound predicts.

On one desktop machine, with -O2, the C++ program needs about \(0.29\) s for the target, the Java port about \(0.34\) s and the Python version about \(4\) s; the three C++ checkpoints add \(2.3\) ms. Of the C++ time the four loops take about \(0.26\) s and the sieve under \(0.01\) s. The access patterns make that split plausible: the sieve walks its table in arithmetic progressions, while each quadruple reads it at an unpredictable position. The seconds themselves depend on the machine.

The method was checked well beyond those checkpoints. A brute force over all positive integer matrices, grouping the squares and keeping those reached twice, agrees with the formula for every \(2\le N\le400\) and at \(N=1000\), \(1500\) and \(2000\). At \(N=1000\) the \(1019\) matrices rebuilt from the admissible quadruples are, as a set, exactly the \(1019\) found by brute force, which tests the whole construction of Section 10 at that bound, not only the total. No matrix with trace below \(1000\) has three positive integer square roots, as Section 3 predicts; in every counted pair \(\det A\) is a perfect square, the two traces satisfy \(t^2=S\pm2\sqrt{\det A}\), and the two roots have opposite determinants \(\mp(p^2-q^2)n^2/4\). The three programs agree with one another at \(F(10^4)=25\,989\), \(F(10^5)=518\,497\) and \(F(10^6)=9\,083\,306\) and at the requested value.

One limitation of the C++ brute force is worth recording, since it is a check rather than part of the method: its 64-bit key holds each matrix entry in sixteen bits, and the entry \(b(a+d)\) first exceeds \(65\,535\) at the bound \(3873\), from where on different matrices could collide in the map. It is only ever called with \(200\), where the margin is wide, but it must not be reused for larger bounds unchanged.

Footnotes and References

The references below supply the background facts used above. The quadratic for \(t^2\), the parameterisation by \((p,q,n,m,u,v)\), the bijection of Section 10 and the operation counts are derived explicitly in this article.

1. Project Euler 420 — 2×2 Positive Integer Matrix. The official statement defines a positive integer matrix as one whose elements are all positive integers, exhibits a matrix that is a square in two different ways, defines F(N) as the number of 2×2 positive integer matrices with trace less than N that can be expressed as a square of a positive integer matrix in two different ways, states F(50) = 7 and F(1000) = 1019, and asks for F(10⁷). The two checkpoints test an implementation without disclosing the requested value.
2. Wikipedia — Cayley–Hamilton theorem. For a 2×2 matrix the page gives the characteristic polynomial p(λ) = λ² − (a + d)λ + (ad − bc) and the resulting identity A² − (a + d)A + (ad − bc)I = 0, which is the form of the square used in Section 1.
3. Wikipedia — Coprime integers. The page lists the property used twice in Sections 5 and 6: if a and b are coprime and a divides the product bc, then a divides c, a generalisation of Euclid's lemma. It also records that the numerator and denominator of a reduced fraction are coprime, which is how p and q are obtained from the two traces.
4. Wikipedia — Divisor function. The page defines σ₀(n), also written τ(n) or d(n), as the number of positive divisors of n, and gives the product formula over the exponents of the prime factorisation. Section 10 uses τ to count the ordered factorisations uv = ℓ.
5. OEIS A000005 — number of divisors of n. The sequence begins 1, 2, 2, 3, 2, 4, 2, 4, 3, 4; in particular τ(1) = 1, τ(4) = 3, τ(8) = 4 and τ(9) = 3, the values used in the worked examples of Sections 12 and 13.
6. cppreference — std::gcd. Declared in <numeric> since C++17, it returns the greatest common divisor of the absolute values of its two integer arguments, and zero only when both are zero. The C++ program uses it to keep the coprime pairs (p, q).
7. cppreference — Fundamental types. The page states that two's complement is the only signed integer representation allowed by the standard as of C++20, with the range −2^(N−1) to 2^(N−1) − 1. This is what makes the bitwise parity test on a negative m well defined.
8. Python documentation — math. The number-theoretic section lists isqrt(n), the integer square root of a non-negative integer n, and gcd(*integers), the greatest common divisor of the arguments; the Python solution uses both.
9. The Java Language Specification, §4.2 — Primitive Types and Values. int and long are 32-bit and 64-bit signed two's-complement integers, with ranges −2147483648 to 2147483647 and −9223372036854775808 to 9223372036854775807. The two's-complement rule makes the Java parity test behave like the C++ one, and the long accumulator holds the growing sum.
10. Wikipedia — Harmonic number. The page gives the asymptotic expansion Hₙ ~ ln n + γ + 1/(2n) − 1/(12n²) + ⋯ with γ ≈ 0.5772156649, which turns the sieve's divisor count into L ln L + O(L).