Detailed mathematical approach

Problem Summary

We must count positive-integer matrices

$A \in M_2(\mathbb{Z}_{>0})$

with

$\operatorname{tr}(A) \lt N$

such that \(A\) has at least two distinct positive-integer square roots. In other words, there exist two different matrices \(B_1\) and \(B_2\) with positive integer entries for which

$B_1^2 = A = B_2^2.$

A direct search over all possible matrices is far too expensive when \(N = 10^7\), so the solution reorganizes the problem around the arithmetic structure shared by two different roots of the same matrix.

Mathematical Approach

Step 1: Start from Two Distinct Square Roots

Let

$B_1=\begin{pmatrix}a & b\\ c & d\end{pmatrix},\qquad B_2=\begin{pmatrix}e & f\\ g & h\end{pmatrix},$

where all eight variables are positive integers and \(B_1 \ne B_2\). Squaring both matrices gives

$B_1^2=\begin{pmatrix}a^2+bc & b(a+d)\\ c(a+d) & d^2+bc\end{pmatrix},$

$B_2^2=\begin{pmatrix}e^2+fg & f(e+h)\\ g(e+h) & h^2+fg\end{pmatrix}.$

Because both squares are the same matrix \(A\), their corresponding entries must agree. Introduce the trace-like and difference-like quantities

$t_1=a+d,\qquad t_2=e+h,\qquad r_1=a-d,\qquad r_2=e-h.$

Then the off-diagonal entries satisfy

$b\,t_1=f\,t_2=A_{12},\qquad c\,t_1=g\,t_2=A_{21},$

and subtracting the two diagonal entries in each square gives

$a^2-d^2=e^2-h^2 \iff r_1 t_1=r_2 t_2.$

If \(t_1=t_2\), then the off-diagonal equalities force \(b=f\) and \(c=g\), and the relation \(r_1 t_1=r_2 t_2\) forces \(r_1=r_2\), hence \(a=e\) and \(d=h\). So distinct roots must have different trace sums. We may therefore assume

$t_1 \lt t_2.$

Step 2: Reduce the Ratio of the Two Trace Sums

Write the ratio \(t_1:t_2\) in lowest terms:

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

Now use \(r_1 t_1=r_2 t_2\). Since \(q r_1 = p r_2\) and \(p,q\) are coprime, there exists an integer \(m\) such that

$r_1=pm,\qquad r_2=qm.$

Therefore the diagonal entries of the two roots are completely determined by \(p,q,n,m\):

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

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

The off-diagonal equalities show that both \(A_{12}\) and \(A_{21}\) are divisible by both \(qn\) and \(pn\). Since \(\gcd(p,q)=1\), they are divisible by \(pqn\). So we may write

$A_{12}=pqn\,u,\qquad A_{21}=pqn\,v,\qquad u,v\in\mathbb{Z}_{>0}.$

Substituting back gives the individual off-diagonal entries of the two roots:

$b=pu,\qquad c=pv,\qquad f=qu,\qquad g=qv.$

Step 3: Positivity and Parity Constraints

The entries \(a,d,e,h\) must be positive integers. Positivity gives

$qn \gt p|m|,$

or equivalently

$|m| \lt \frac{qn}{p}.$

Integrality imposes parity conditions, because each diagonal entry is a half-integer expression. We need \(qn \pm pm\) and \(pn \pm qm\) to be even. Since \(p\) and \(q\) are coprime, this reduces to two cases:

If \(p\) and \(q\) have opposite parity, then both \(n\) and \(m\) must be even.

If \(p\) and \(q\) are both odd, then \(m \equiv n \pmod{2}\).

These are exactly the parity filters used by the implementations when they enumerate valid values of \(m\).

Step 4: Determine the Product \(uv\)

Set

$\ell=uv.$

Now compare the upper-left entries of \(B_1^2\) and \(B_2^2\):

$\frac{(qn+pm)^2}{4}+p^2\ell=\frac{(pn+qm)^2}{4}+q^2\ell.$

After expanding and simplifying, the mixed terms cancel and we obtain

$4\ell=n^2-m^2.$

Hence

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

Because \(u\) and \(v\) are positive, we need \(\ell \gt 0\), so automatically \(n^2 \gt m^2\). The same identity also makes the lower-right entries agree, so this single equation is the key Diophantine compatibility condition.

Step 5: Why the Divisor Function Appears

For fixed \(p,q,n,m\), the diagonal parts of both roots are fixed, and so is the common scale factor \(pqn\) in the off-diagonal entries of \(A\). The only remaining freedom is to factor \(\ell\) as

$\ell=uv,\qquad u,v\in\mathbb{Z}_{>0}.$

Each positive divisor \(u\mid \ell\) determines exactly one ordered pair

$\left(u,\frac{\ell}{u}\right),$

and therefore exactly one matrix \(A\). The number of such ordered factor pairs is the divisor function \(\tau(\ell)\). That is why every admissible quadruple \((p,q,n,m)\) contributes

$\tau\!\left(\frac{n^2-m^2}{4}\right)$

to the final count.

Step 6: Trace Bound and Final Summation Formula

Using \(B_1\), we can compute the trace of \(A\):

$\operatorname{tr}(A)=a^2+d^2+2bc.$

Substitute \(a,d,b,c\) from the parameterization:

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

Now replace \(\ell\) by \(\frac{n^2-m^2}{4}\):

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

So the trace condition \(\operatorname{tr}(A)\lt N\) becomes

$\frac{(p^2+q^2)n^2}{2} \lt N,$

or

$n \le \sqrt{\frac{2N-1}{p^2+q^2}}.$

This gives the final counting formula:

$F(N)=\sum_{\substack{p \gt q \ge 1\\ \gcd(p,q)=1}} \ \sum_{n=1}^{\left\lfloor\sqrt{\frac{2N-1}{p^2+q^2}}\right\rfloor} \ \sum_{m}\tau\!\left(\frac{n^2-m^2}{4}\right),$

where \(m\) runs over the integers satisfying

$|m| \lt \frac{qn}{p}$

together with the parity conditions from the previous step. This is the exact formula used by the implementations.

Step 7: Small Worked Example

Take

$p=2,\qquad q=1,\qquad n=4,\qquad m=0.$

Then

$\ell=\frac{4^2-0^2}{4}=4,\qquad \operatorname{tr}(A)=\frac{(2^2+1^2)\cdot 4^2}{2}=40.$

Since \(\tau(4)=3\), this single admissible quadruple contributes three matrices, corresponding to the factor pairs \((u,v)=(1,4),(2,2),(4,1)\). For \((u,v)=(1,4)\), the two roots are

$\begin{pmatrix}2 & 2\\ 8 & 2\end{pmatrix},\qquad \begin{pmatrix}4 & 1\\ 4 & 4\end{pmatrix},$

and both square to

$\begin{pmatrix}20 & 8\\ 32 & 20\end{pmatrix}.$

This example shows exactly how one admissible arithmetic state can generate several matrices through the divisors of \(\ell\).

How the Code Works

The C++, Python, and Java implementations follow the same structure. They first derive an upper bound for \(\ell=\frac{n^2-m^2}{4}\). The smallest possible value of \(p^2+q^2\) is \(5\), attained at \((p,q)=(2,1)\), so \(n^2 \le \frac{2N-1}{5}\), and therefore \(\ell \le \frac{n^2}{4}\). With that bound in hand, the implementation precomputes \(\tau(x)\) for every needed \(x\) by a divisor sieve.

After the sieve, the program iterates over coprime pairs \((p,q)\) with \(p \gt q\), computes the admissible range of \(n\) from the trace bound, and then scans the valid values of \(m\) respecting the positivity and parity conditions. For each admissible state it evaluates \(\ell=(n^2-m^2)/4\) and adds \(\tau(\ell)\) to the total. The implementations also verify the reduction with small checkpoints; in particular, they confirm \(F(50)=7\) and \(F(1000)=1019\), and the C++ version cross-checks the formula against direct enumeration for a small bound.

Complexity Analysis

The divisor sieve up to

$\left\lfloor \frac{1}{4}\left\lfloor \sqrt{\frac{2N-1}{5}} \right\rfloor^2 \right\rfloor$

uses \(O(L\log L)\) time and \(O(L)\) memory, where \(L\) denotes that bound. Since \(L\) is proportional to \(N\), this is effectively \(O(N\log N)\) time and \(O(N)\) memory for the preprocessing stage.

The arithmetic enumeration over \((p,q,n,m)\) is far smaller than a brute-force search over all positive roots. Its aggregate size is also \(O(N\log N)\), so the whole method remains practical for \(N=10^7\). The crucial gain is that the program never enumerates matrices \(A\) directly; it only enumerates reduced number-theoretic states that are known to produce at least two square roots.

References

1. Problem page: https://projecteuler.net/problem=420
2. Divisor function: Wikipedia — Divisor function
3. Matrix square root: Wikipedia — Square root of a matrix
4. Coprime integers: Wikipedia — Coprime integers
5. Greatest common divisor: Wikipedia — Greatest common divisor