Detailed mathematical approach

Problem Summary

We must count integer pairs \((a,b)\neq(0,0)\) such that

$a^2+ab+41b^2\le N.$

This binary quadratic form is positive definite, so all valid pairs lie in a bounded region. A naive two-dimensional scan would still be far too slow for \(N=10^{16}\). The key observation is that a simple linear change of variables turns the form into an axis-aligned ellipse with a parity condition, which makes the counting almost one-dimensional.

Mathematical Approach

Let \(F(N)\) denote the number of nonzero integer pairs \((a,b)\) satisfying \(a^2+ab+41b^2\le N\). The implementations count \(F(N)\) by rewriting the quadratic form and then sweeping over one coordinate.

Step 1: Diagonalize the quadratic form

Introduce the substitution

$x=2a+b,\qquad y=b.$

Then

$4(a^2+ab+41b^2)=(2a+b)^2+163b^2=x^2+163y^2.$

So the original inequality is equivalent to

$x^2+163y^2\le 4N.$

The coefficient \(163\) is the absolute value of the discriminant \(1-4\cdot 41=-163\), so its appearance is natural rather than accidental.

Step 2: Recover the integrality condition

The substitution is reversible:

$a=\frac{x-y}{2},\qquad b=y.$

Therefore \((a,b)\) is an integer pair exactly when \(x-y\) is even, which is the same as

$x\equiv y\pmod{2}.$

This gives the equivalent counting formula

$F(N)=\#\left\{(x,y)\in\mathbb{Z}^2:x^2+163y^2\le 4N,\ x\equiv y\pmod{2}\right\}-1.$

The final \(-1\) removes the origin, which is the unique lattice point producing the value \(0\).

Step 3: Bound the outer variable

For any admissible point,

$163y^2\le 4N.$

Hence

$|y|\le Y=\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor.$

So there are only \(2Y+1\) possible values of \(y\). This collapses the search from a two-variable problem to a one-variable sweep over horizontal slices of the ellipse.

Step 4: Count the admissible \(x\)-values for fixed \(y\)

Fix one integer \(y\) with \(|y|\le Y\), and define

$R(y)=4N-163y^2,\qquad M(y)=\left\lfloor\sqrt{R(y)}\right\rfloor.$

Then \(x^2\le R(y)\) is equivalent to

$-M(y)\le x\le M(y).$

Without the parity restriction, this interval contains \(2M(y)+1\) integers. We only want those whose parity matches \(y\).

Step 5: Convert parity into a closed formula

The interval \([-M,M]\) is symmetric around \(0\). If \(M\) is even, it contains \(M+1\) even integers and \(M\) odd integers. If \(M\) is odd, the roles reverse. Therefore the number of valid \(x\)-values for a fixed \(y\) is

$C(y)=\begin{cases} M(y)+1,& y\equiv M(y)\pmod{2},\\ M(y),& y\not\equiv M(y)\pmod{2}. \end{cases}$

This is the crucial simplification: one integer square root and one parity test replace a full inner enumeration of \(x\).

Step 6: Sum all slices

Adding the contribution of every permitted \(y\) gives

$F(N)=\sum_{y=-Y}^{Y} C(y)-1.$

Because \(M(y)=M(-y)\), the picture is symmetric. The implementations still loop over the full range directly, which keeps the logic simple and avoids any special handling at \(y=0\).

Worked Example: \(N=1000\)

For \(N=1000\),

$Y=\left\lfloor\sqrt{\frac{4000}{163}}\right\rfloor=4.$

Now evaluate the slices:

$\begin{aligned} y=0 &:\quad R=4000,\quad M=63,\quad C=63,\\ y=\pm 1 &:\quad R=3837,\quad M=61,\quad C=62,\\ y=\pm 2 &:\quad R=3348,\quad M=57,\quad C=57,\\ y=\pm 3 &:\quad R=2533,\quad M=50,\quad C=50,\\ y=\pm 4 &:\quad R=1392,\quad M=37,\quad C=37. \end{aligned}$

Therefore

$F(1000)=63+2(62+57+50+37)-1=474.$

This matches the count returned by the implementation and confirms the transformation and parity formula.

How the Code Works

The C++, Python, and Java implementations follow the same mathematical plan. First they compute an exact floor square root. Python uses the standard exact integer routine, while the C++ and Java versions start from a floating approximation and then correct it upward or downward until the exact integer floor is reached. That correction step matters because the input is large enough that a one-unit error near a square boundary would change the final count.

After that, the implementation finds the largest permitted absolute value of the outer coordinate, iterates through every integer in that range, forms the remaining budget \(4N-163y^2\), takes its integer square root, and applies the parity formula above to count the matching values of the inner coordinate. Each loop iteration contributes one horizontal slice of the ellipse, and the origin is removed once at the end.

Complexity Analysis

The loop length is

$2\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor+1=\Theta(\sqrt N).$

Each iteration performs constant-time arithmetic plus one integer square root, so the total running time is \(O(\sqrt N)\). Only a handful of integer variables are stored, so the memory usage is \(O(1)\).

Footnotes and References

1. Problem page: Project Euler 804
2. Binary quadratic forms: Wikipedia - Binary quadratic form
3. Discriminant: Wikipedia - Discriminant
4. Integer square root: Wikipedia - Integer square root
5. Lattice points: Wikipedia - Lattice point