Detailed mathematical approach

Problem Summary

Write

$r=AO=CO,\qquad s=BO=DO=\frac{BD}{2},\qquad 0<r\le s.$

Then every biclinic quadrilateral lives on one fixed pair of concentric circles centered at \(O\): the vertices \(A,C\) lie on radius \(r\), while \(B,D\) are the opposite endpoints of a diameter of radius \(s\).

The key invariant is

$AB^2+BC^2+CD^2+AD^2=4(r^2+s^2).$

So if we define

$n=r^2+s^2,$

then the condition in the problem is simply \(n\le N/4\). The whole problem becomes: for each \(n\), how many biclinic quadrilaterals are produced by the representations of \(n\) as a sum of two squares?

Mathematical Approach

1) Every side pair gives a strict representation of the same \(n\).

Look at vertex \(A\). Let

$u=AB,\qquad v=AD,\qquad u<v.$

Define

$x=\frac{u+v}{2},\qquad y=\frac{v-u}{2},\qquad u=x-y,\quad v=x+y.$

In triangle \(BAD\), the point \(O\) is the midpoint of \(BD\). Apollonius gives

$u^2+v^2=2(AO^2+BO^2)=2(r^2+s^2)=2n.$

But also

$u^2+v^2=(x-y)^2+(x+y)^2=2(x^2+y^2).$

Therefore

$x^2+y^2=n.$

Because \(u\) and \(v\) are distinct positive integers, we get a strict representation \(x>y>0\). The same argument at vertex \(C\) gives another strict representation of the same \(n\).

2) The reverse construction.

Suppose we already know one representation

$n=s^2+r^2,\qquad s\ge r>0,$

which will be used for the diagonal data \(BD=2s\) and \(AO=CO=r\). Now take another strict representation

$n=x^2+y^2,\qquad x>y>0.$

Set

$u=x-y,\qquad v=x+y.$

Then

$u^2+v^2=2n=2(r^2+s^2),$

so by Apollonius any point \(P\) with \(PB=u\) and \(PD=v\) automatically satisfies \(PO=r\). Thus the pair \((x,y)\) determines a valid vertex on the inner circle.

To make the circles intersect, we need the triangle inequalities with base \(BD=2s\):

$v-u=2y<2s<2x=u+v.$

This is exactly why the diagonal representation must be the smallest one in the ordering below.

3) Order the representations of \(n\).

Let

$\mathcal R(n)=\{(x,y)\in\mathbb Z_{>0}^2:\ x>y,\ x^2+y^2=n\},\qquad m(n)=|\mathcal R(n)|.$

Sort the pairs by increasing \(x\):

$x_1<x_2<\cdots<x_{m(n)}.$

Because \(x_i^2+y_i^2=n\), the corresponding \(y_i\) values strictly decrease:

$y_1>y_2>\cdots>y_{m(n)}.$

On the branch \(x>\sqrt{n/2}\), the two functions

$x-y\quad\text{and}\quad x+y$

behave monotonically in opposite directions:

$x-y\ \text{increases},\qquad x+y\ \text{decreases}.$

So if we choose three strict representations

$ (x_i,y_i),\ (x_j,y_j),\ (x_k,y_k)\qquad (i<j<k),$

then the smallest one \((x_i,y_i)\) is the only possible diagonal pair \((s,r)\), because it guarantees \(x_j>s\) and \(x_k>s\), while the reverse choice would violate the triangle inequality.

4) One 3-subset of representations gives one quadrilateral.

Take the diagonal data from the smallest representation:

$s=x_i,\qquad r=y_i,\qquad BD=2x_i,\qquad AO=CO=y_i.$

Use the next two representations for the two vertices:

$AB=x_j-y_j,\qquad AD=x_j+y_j,$

$BC=x_k-y_k,\qquad CD=x_k+y_k.$

The monotonicity above gives the strict side order automatically:

$AB<BC<CD<AD.$

Hence every 3-element subset of \(\mathcal R(n)\) contributes exactly one biclinic integral quadrilateral. This yields the generic term

$\binom{m(n)}{3}.$

5) The special case \(n=2t^2\).

If

$n=2t^2,$

then besides the strict representations there is also the symmetric one \((t,t)\). It is not counted in \(m(n)\) because it is not strict, but it is perfectly valid for the diagonal data:

$s=t,\qquad r=t,\qquad AO=BO=CO=DO=t.$

Once this diagonal pair is fixed, any two strict representations of \(n\) can play the roles of the two vertices, so this case contributes

$\binom{m(n)}{2}.$

6) Final counting formula.

Putting the generic and special cases together,

$B(N)=\sum_{1\le n\le N/4}\left[\binom{m(n)}{3}+\mathbf 1_{\,n=2t^2}\binom{m(n)}{2}\right].$

Worked Examples

Official sample. The example in the statement has

$AO=CO=23,\qquad BO=DO=24,$

so

$n=23^2+24^2=1105.$

The strict representations of \(1105\) are

$1105=24^2+23^2=31^2+12^2=32^2+9^2=33^2+4^2.$

Choose \((24,23)\) as the diagonal pair. Then using \((31,12)\) and \((33,4)\) for the two vertices gives

$AB=31-12=19,\qquad AD=31+12=43,$

$BC=33-4=29,\qquad CD=33+4=37,$

and

$BD=2\cdot 24=48,\qquad AO=CO=23.$

This is exactly the quadrilateral from the problem statement. Since \(m(1105)=4\), the value \(n=1105\) actually contributes \(\binom{4}{3}=4\) quadrilaterals in total.

Why the extra \(\binom{m(n)}{2}\) term exists. Consider

$1250=25^2+25^2=31^2+17^2=35^2+5^2.$

Here \(m(1250)=2\) because only the two strict representations are counted. The generic term is \(\binom{2}{3}=0\), but the symmetric pair \((25,25)\) can be used as the diagonal representation, so we get one extra quadrilateral:

$\binom{2}{2}=1.$

Algorithm

1) Work blockwise in \(n\). Let \(T=N/4\). The code processes intervals

$[L,H]\subseteq [1,T]$

of width \(W\). For one block it stores a local counter array `counts[n-L]`.

2) Enumerate all strict representations in the block.

For every \(y\ge 1\), the smallest possible strict value is

$y^2+(y+1)^2.$

If that already exceeds \(H\), we stop. Otherwise the feasible \(x\)-range is

$x_{\min}=\max\!\left(y+1,\left\lceil \sqrt{L-y^2}\right\rceil\right),\qquad x_{\max}=\left\lfloor \sqrt{H-y^2}\right\rfloor.$

Every pair \((x,y)\) with \(x_{\min}\le x\le x_{\max}\) contributes one strict representation to

$n=x^2+y^2.$

3) Convert representation counts into quadrilateral counts.

After the block is filled, every touched value \(n\) has \(m(n)\) stored in its counter. The code adds

$\binom{m(n)}{3}$

and, if \(n/2\) is a square, also

$\binom{m(n)}{2}.$

Then the touched counters are reset to zero and the next block is processed.

4) Why the implementation is fast.

It never factors every integer \(n\). It only enumerates lattice points \((x,y)\) that actually satisfy \(x^2+y^2\le T\). The `touched` list ensures that only entries hit in the current block are visited when the block is finalized.

Complexity Analysis

Let \(T=N/4\). The memory usage is \(O(W)\) for block width \(W\). The running time is essentially proportional to the number of strict lattice representations

$x^2+y^2\le T,\qquad x>y>0,$

that are generated across all blocks. The code also parallelizes naturally because different blocks are independent.

Checks And Final Result

The C++ implementation verifies

$B(10\,000)=49,\qquad B(1\,000\,000)=38239,$

and the default full run returns

$B(10\,000\,000\,000)=2466018557.$

Further Reading

1. Problem page: https://projecteuler.net/problem=311
2. Apollonius's theorem: https://en.wikipedia.org/wiki/Apollonius's_theorem
3. Sum of two squares theorem: https://en.wikipedia.org/wiki/Fermat's_theorem_on_sums_of_two_squares