Detailed mathematical approach

Problem Summary

The solution works with a canonical parametrization by primitive integer pairs \((p,q)\) with \(q\ge 0\) and \(\gcd(p,q)=1\). In that parametrization, the size attached to a solution is the discriminant-\(13\) binary quadratic form

$z=\left|3p^2-7pq+3q^2\right|,$

and we must count all admissible primitive pairs for which \(z\le N\). A second quadratic coordinate must stay positive as well, so the ratio \(p/q\) is not free: only three regions survive after the sign and symmetry reductions.

The isolated case \(q=0\) gives \(z=3\). Every other contribution comes from \(q\ge 1\), which is why the implementations loop over \(q\) and count valid \(k\)-intervals instead of scanning the original search space directly.

Mathematical Approach

Introduce the two quadratic forms

$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$

The counted size is \(|Q(p,q)|\), while admissibility requires \(X(p,q)>0\), \(|p|>q\), and \(\gcd(p,q)=1\). The code turns those conditions into interval counts in a new variable \(k\).

Step 1: Replace \(p\) by a positive offset \(k\)

Because \(|p|>q\), every admissible pair with \(q>0\) can be written in exactly one of the forms

$p=q+k,\qquad p=-q-k.$

with \(k\ge 1\). Coprimality is preserved under this substitution:

$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$

So after fixing \(q\), the arithmetic condition becomes simply \(\gcd(q,k)=1\).

Step 2: Split the search into the three surviving regions

For \(p=q+k\), the auxiliary form becomes

$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q).$

If we write

$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3,$

then \(X>0\) forces either \(k<\alpha q\) or \(k>\beta q\). These are the two branches called region A and region D:

$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$

$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$

For \(p=-q-k\), we get

$X(-q-k,q)=13q^2+8qk+k^2>0,$

so positivity is automatic, and the size becomes

$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2.$

Thus the entire count is the sum of region A, region C, region D, and the single base solution \(z=3\).

Step 3: Count coprime \(k\) by Möbius inversion

For fixed \(q\), define

$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}.$

Using inclusion-exclusion over the prime divisors of \(q\),

$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor.$

Only squarefree divisors matter, so after factoring \(q\) once, the implementations generate all \(2^{\omega(q)}\) squarefree divisors together with their Möbius signs. Any interval \([L,U]\) is then counted by

$C_q(U)-C_q(L-1).$

Step 4: Remove the non-canonical branch modulo \(13\)

The discriminant-\(13\) form has a simple congruence:

$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$

Therefore

$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$

The implementations keep only the canonical branch with \(13\nmid Q(p,q)\), so one residue class must be removed:

$p=q+k \Rightarrow k\equiv -2q \pmod{13},$

$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$

When \(q\equiv 0\pmod{13}\), coprimality already forbids \(k\equiv 0\pmod{13}\), so there is no extra subtraction. Otherwise the forbidden class is counted with the same Möbius table, but now on an arithmetic progression modulo \(13\).

Step 5: Turn the size bound \(z\le N\) into explicit endpoints

Each region contributes an interval in \(k\).

Region A starts with \(1\le k\le \lfloor \alpha q\rfloor\). The polynomial \(z_A(q,k)=q^2+qk-3k^2\) is concave, so the inequality \(z_A\le N\) can fail only on a middle interval. Solving \(z_A=N\) gives

$k=\frac{q\pm\sqrt{13q^2-12N}}{6}.$

If \(13q^2\le 12N\), the whole region survives; otherwise the integers between the two roots must be removed.

Region C is monotone increasing, so its upper bound comes from

$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}.$

Region D is also monotone once \(k>\beta q\), so it uses

$k\ge \lfloor \beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}.$

The implementations compute these bounds from integer square roots and then adjust the endpoints by a few exact checks, eliminating rounding errors.

Worked Example: \(N=100\) and \(q=1\)

Here \(\lfloor \alpha\rfloor=0\), so region A is empty.

Region C satisfies

$13+13k+3k^2\le 100,$

which gives \(k=1,2,3\), hence the contributions \(29,51,79\).

Region D starts at \(k=\lfloor \beta\rfloor+1=4\). There we need

$3k^2-k-1\le 100,$

so \(k=4,5\), giving \(43\) and \(69\).

For \(q=1\), the forbidden classes are \(k\equiv 11\pmod{13}\) in regions A and D, and \(k\equiv 0\pmod{13}\) in region C, so no subtraction occurs in this layer. Together with the isolated case \(q=0\), this shows exactly how the counting proceeds.

How the Code Works

The C++, Python, and Java implementations first set \(Q_{\max}=\lfloor\sqrt N\rfloor\) and build a smallest-prime-factor table up to that limit. This lets them factor every \(q\) quickly and generate all squarefree divisors needed for the Möbius sums.

For each \(q\), the implementation evaluates the three regions separately. Region A is handled as a prefix count up to \(\lfloor\alpha q\rfloor\) minus the middle interval where \(z_A>N\). Regions C and D are monotone, so they are counted with direct prefix differences. In every case the coprime filter and, when necessary, the forbidden residue class modulo \(13\) are applied through the same divisor table.

The C++ implementation additionally splits the outer \(q\)-range across worker threads and accumulates partial sums in parallel. The Python and Java implementations perform the same arithmetic serially. After all \(q\ge 1\) contributions are finished, the isolated primitive case \(z=3\) is added once \(N\ge 3\).

Complexity Analysis

Let \(Q=\lfloor\sqrt N\rfloor\). Building the smallest-prime-factor sieve costs \(O(Q\log\log Q)\) time and \(O(Q)\) memory. For a fixed \(q\), the number of squarefree divisors is \(2^{\omega(q)}\), and each region uses only a constant number of sums over that table.

So the total running time is

$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right),$

which is close to \(O(Q\log Q)\) on average and behaves near-linearly in practice. Memory usage stays \(O(Q)\), with only small per-\(q\) temporary arrays besides the sieve.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=769
2. Binary quadratic form: Wikipedia — Binary quadratic form
3. Möbius inversion formula: Wikipedia — Möbius inversion formula
4. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
5. Greatest common divisor: Wikipedia — Greatest common divisor