Solution

Answer: 2611227421428

For a quadrilateral with side lengths $a,b,c,d$, the maximal possible area is attained by a cyclic quadrilateral. By Brahmagupta’s formula,

$$K^2=(s-a)(s-b)(s-c)(s-d), \qquad s=\frac{a+b+c+d}{2}.$$

So the problem reduces to studying integer solutions of this expression.

Define

$$U=s-a,\quad V=s-b,\quad W=s-c,\quad T=s-d.$$

Then

$$U\ge V\ge W\ge T>0$$

because $a\le b\le c\le d$, and

$$a=s-U,\ b=s-V,\ c=s-W,\ d=s-T.$$

Also,

$$U+V+W+T=2s=a+b+c+d.$$

The maximal area becomes

$K^2=UVWT$

Hence $K$ is an integer iff $UVWT$ is a perfect square.

The side lengths are integers iff $U,V,W,T$ all have the same parity, equivalently

$$U+V+W+T \equiv 0 \pmod 2.$$

Finally, $a>0$ gives

$$s-U>0 \iff U<V+W+T.$$

So we must enumerate all ordered quadruples

$$U\ge V\ge W\ge T>0$$

such that

$UVWT=k^2$ for some $1\le k\le n$,
$U<V+W+T$,
$U+V+W+T$ is even.

For every valid quadruple, the contribution to $SP(n)$ is

$$U+V+W+T.$$

Efficient enumeration

Fix $k\le n$. We enumerate divisors of $k^2$.

Let

$$UVWT=k^2.$$

Choose $T\mid k^2$, then $W\mid k^2/T$, and define

$$R=\frac{k^2}{TW}=UV.$$

Now $V$ must divide $R$, and

$$U=\frac{R}{V}.$$

The ordering condition $U\ge V$ implies

$$V\le \sqrt R.$$

The positivity condition becomes

$$\frac{R}{V}<V+W+T.$$

This gives a lower bound for $V$, dramatically pruning the search.

The resulting algorithm runs comfortably fast for $n=10^6$.

Python implementation

from math import isqrt
from bisect import bisect_left

def build_spf(limit):
    """Smallest prime factor sieve."""
    spf = list(range(limit + 1))

    for i in range(2, isqrt(limit) + 1):
        if spf[i] == i:
            for j in range(i * i, limit + 1, i):
                if spf[j] == j:
                    spf[j] = i

    return spf

def factorize(n, spf):
    """Prime factorization of n."""
    out = []

    while n > 1:
        p = spf[n]
        e = 0

        while n % p == 0:
            n //= p
            e += 1

        out.append((p, e))

    return out

def divisors_of_k2_upto_k(k, fac):
    """
    Generate all divisors d of k^2 with d <= k.
    """

    divs = [1]

    for p, e in fac:
        current = divs[:]

        pe = 1
        for _ in range(2 * e):
            pe *= p

            for d in current:
                v = d * pe
                if v <= k:
                    divs.append(v)

    divs.sort()
    return divs

def SP(n):
    spf = build_spf(n)

    total = 0

    for k in range(1, n + 1):

        k2 = k * k

        fac = factorize(k, spf)

        divs = divisors_of_k2_upto_k(k, fac)

        dlen = len(divs)

        for ti, T in enumerate(divs):

            k2_div_T = k2 // T

            for wi in range(ti, dlen):

                W = divs[wi]

                if W > k2_div_T:
                    break

                if k2_div_T % W:
                    continue

                R = k2_div_T // W   # R = U*V

                vmax = isqrt(R)

                if vmax < W:
                    continue

                S = W + T

                # Need:
                # R/V < V + S
                #
                # Solve:
                # V^2 + S*V - R > 0

                disc = S * S + 4 * R
                root = isqrt(disc)

                vmin = (root - S) // 2 + 1

                if vmin < W:
                    vmin = W

                if vmin > vmax:
                    continue

                lo = bisect_left(divs, vmin)

                for vi in range(lo, dlen):

                    V = divs[vi]

                    if V > vmax:
                        break

                    if R % V:
                        continue

                    U = R // V

                    # positivity condition
                    if U >= V + S:
                        continue

                    perimeter = U + V + S

                    # all U,V,W,T same parity
                    if perimeter % 2 == 0:
                        total += perimeter

    return total

# verification
assert SP(10) == 186
assert SP(100) == 23238

print(SP(1_000_000))

Code walkthrough

1. Convert the geometry problem to arithmetic

Using Brahmagupta’s formula, we replace the sides by

$$(U,V,W,T)=(s-a,s-b,s-c,s-d).$$

Then:

area condition becomes $UVWT=k^2$,
ordering becomes $U\ge V\ge W\ge T$,
side positivity becomes $U<V+W+T$,
integrality becomes even perimeter.

2. Enumerate divisors

For each $k\le n$, we factor $k$ using a smallest-prime-factor sieve.

From the factorization of $k$, we generate all divisors of $k^2$ up to $k$.

3. Construct $(U,V,W,T)$

We choose:

$T$,
$W$,
then $V$,

and compute

$$U=\frac{k^2}{VWT}.$$

Bounds derived from the inequality

$$U<V+W+T$$

greatly reduce the search space.

4. Verify examples

The code checks:

assert SP(10) == 186
assert SP(100) == 23238

matching the statement exactly.

Final verification

The computation gives

$$SP(1,000,000)=2611227421428.$$

The value is positive and of the expected magnitude (roughly quadratic growth compared with the $n=100$ case), and all provided sample values are reproduced exactly.

Answer: 2611227421428

Project Euler Problem 681