Detailed mathematical approach

Problem Summary

Define \(R(k)\) to be the number of representations of \(k\) in the form

$k=a^2+3ab+b^2,\qquad a \gt b \gt 0.$

The problem asks for

$f(n,r)=\#\{k\le n : R(k)=r\},$

and in particular for \(f(10^{15},40)\). A brute-force search over pairs \((a,b)\) is hopeless, so the solution classifies represented integers by their prime-factor structure in the quadratic field of discriminant \(5\).

Mathematical Approach

The key observation is that the form behaves like a norm. Once that structure is exposed, the exact representation count depends only on the primes that split modulo \(5\), while the factors \(5^t\) and inert squares merely scale the represented integer without changing its multiplicity.

Step 1: Rewrite the Form as a Norm

Let

$\alpha=\frac{3+\sqrt5}{2},\qquad \bar{\alpha}=\frac{3-\sqrt5}{2}.$

Then

$a^2+3ab+b^2=(a+b\alpha)(a+b\bar{\alpha}).$

So the quadratic form is the norm of the algebraic integer \(a+b\alpha\) in the field \(\mathbb{Q}(\sqrt5)\). This is why prime splitting modulo \(5\) controls representability: the discriminant of the form is \(5\), and the arithmetic reduces to how rational primes factor in that field.

Step 2: Characterize the Integers That Can Be Represented

For primes \(p\neq 5\), the splitting behavior is:

$p\equiv 1,4\pmod5 \quad \text{split},\qquad p\equiv 2,3\pmod5 \quad \text{inert}.$

Therefore a represented integer has the shape

$k=5^t\left(\prod_{p\equiv 1,4\!\!\!\pmod5} p^{e_p}\right)\left(\prod_{q\equiv 2,3\!\!\!\pmod5} q^{2f_q}\right).$

The inert primes may appear only with even exponent, while the ramified prime \(5\) may appear with any exponent. It is convenient to separate this as

$k=x\cdot 5^t\cdot u^2,$

where

$x=\prod_{p\equiv 1,4\!\!\!\pmod5} p^{e_p}$

is the split core, and every prime factor of \(u\) satisfies \(q\equiv 2,3\pmod5\). The important point is that the exact number of representations depends only on \(x\).

Step 3: Count Representations from the Split Core

For the split core \(x=\prod p^{e_p}\), each split prime \(p\) contributes \(e_p+1\) choices when the factorization is distributed between the two conjugate algebraic factors. Multiplying over all split primes gives

$\tau(x)=\prod_{p\mid x}(e_p+1),$

which is exactly the ordinary divisor function of \(x\).

After restricting to the chamber \(a \gt b \gt 0\), conjugate factorizations are paired, and the middle symmetric case does not create a second interior solution. Hence

$R(k)=\left\lfloor\frac{\tau(x)}{2}\right\rfloor.$

This explains the condition used by the implementations: \(k\) has exactly \(r\) representations if and only if

$\tau(x)\in\{2r,\,2r+1\}.$

Step 4: Turn Exact Multiplicity into Exponent Patterns

Suppose \(\tau(x)=T\), where \(T\) is either \(2r\) or \(2r+1\). Since

$\tau(x)=\prod_{j=1}^s (e_j+1),$

every multiplicative partition

$T=f_1f_2\cdots f_s,\qquad f_j\ge 2,$

produces an exponent pattern

$e_j=f_j-1.$

For example, when \(r=40\), the target values are \(80\) and \(81\). A factorization \(80=5\cdot 4\cdot 4\) yields exponents \((4,3,3)\), corresponding to split cores of the form

$x=p_1^4p_2^3p_3^3,\qquad p_i\equiv 1,4\pmod5,\quad p_i\ \text{distinct}.$

So the first stage of the algorithm is purely combinatorial: generate all exponent multisets whose \((e_j+1)\)-product is \(2r\) or \(2r+1\).

Step 5: Count the Multiplier Part \(5^t u^2\)

Fix one split core \(x\). Any represented number with that same representation count is

$k=x\cdot 5^t\cdot u^2,$

with \(u\) composed only of inert primes. Thus for

$q=\left\lfloor\frac{n}{x}\right\rfloor,$

the number of admissible multipliers is

$G(q)=\#\{5^t u^2\le q : u \text{ uses only inert primes}\}.$

If we define

$I(y)=\#\{u\le y : \text{every prime factor of }u\text{ is }2\text{ or }3\pmod5\},$

then

$G(q)=\sum_{t\ge 0} I\!\left(\left\lfloor\sqrt{\frac{q}{5^t}}\right\rfloor\right),$

where the sum is finite because \(5^t\le q\). This is exactly the table-driven multiplier count built by the implementations.

Step 6: Enumerate Feasible Split Cores Efficiently

For each exponent pattern, the solver assigns those exponents to distinct split primes \(p\equiv 1,4\pmod5\). Because the bound \(x\le n\) matters, different permutations of the same exponent multiset can lead to different feasible products, so distinct orderings must be explored.

The search is pruned aggressively. Two facts are especially useful:

$\text{largest exponents should be paired with the smallest split primes to get the minimal possible product},$

and

$\text{if even that minimal remaining tail already exceeds }n,\text{ the whole branch can be discarded.}$

After all feasible split cores are generated, the final answer is

$f(n,r)=\sum_{x\in \mathcal{X}_{n,r}} G\!\left(\left\lfloor\frac{n}{x}\right\rfloor\right),$

where \(\mathcal{X}_{n,r}\) is the set of split cores satisfying \(\tau(x)\in\{2r,2r+1\}\) and \(x\le n\).

Worked Example

Take

$209=11\cdot 19.$

Both primes are split because \(11\equiv 1\pmod5\) and \(19\equiv 4\pmod5\). Hence the split core is \(x=209\) and

$\tau(x)=(1+1)(1+1)=4,$

so

$R(209)=\left\lfloor\frac{4}{2}\right\rfloor=2.$

Indeed, the two admissible representations are

$209=13^2+3\cdot 13\cdot 1+1^2=8^2+3\cdot 8\cdot 5+5^2.$

For an odd divisor count, consider

$121=11^2,\qquad \tau(121)=3.$

which gives

$R(121)=\left\lfloor\frac{3}{2}\right\rfloor=1,$

and indeed

$121=7^2+3\cdot 7\cdot 3+3^2.$

This is exactly why the target condition is \(\tau(x)\in\{2r,2r+1\}\) rather than only \(2r\).

How the Code Works

The C++, Python, and Java implementations all follow the same decomposition.

First, they generate every exponent multiset whose \((e_j+1)\)-product equals \(2r\) or \(2r+1\), and they immediately discard any pattern whose smallest possible split core already exceeds \(n\). Next, they estimate how large a split prime could ever be in a feasible core, sieve all split primes up to that bound, and precompute the powers needed for every exponent appearing in the surviving patterns.

Separately, they count inert-only integers up to a square-root bound, turn that prefix information into the multiplier table \(G(q)\), and then perform a depth-first search over increasing split primes. Every completed split core \(x\) contributes \(G(\lfloor n/x\rfloor)\). The C++ implementation also evaluates independent pattern-ordering tasks in parallel, while the Python and Java implementations keep the same mathematical structure with lighter orchestration.

Complexity Analysis

Let \(x_{\min}\) be the smallest feasible split core and let \(M=\lfloor n/x_{\min}\rfloor\). Building the inert-only prefix up to \(\lfloor\sqrt{M}\rfloor\) takes sieve-like time near \(O(\sqrt{M}\log\log M)\), and filling the multiplier table costs \(O(M\log M)\) in the loose bound coming from the powers of \(5\). The remaining cost is the DFS over feasible split cores.

In the worst case that search is combinatorial, but for this problem it stays practical because only exponent patterns arising from \(2r\) and \(2r+1\) are considered, the number of distinct exponents is small, and the minimal-tail pruning removes impossible branches very early. Memory usage is dominated by the multiplier table, the inert-only prefix, and the split-prime power tables.

Footnotes and References

1. Problem page: Project Euler 586
2. Binary quadratic forms: Wikipedia - Binary quadratic form
3. Norms in algebraic number theory: Wikipedia - Norm (mathematics)
4. Quadratic fields: Wikipedia - Quadratic field
5. Quadratic residues and splitting modulo \(5\): Wikipedia - Quadratic residue
6. Divisor function: Wikipedia - Divisor function