Solution

Answer: 22472871503401097

Let the four chosen side lengths be $a \le b \le c \le d$.

For a fixed set of side lengths, the quadrilateral with maximal area is cyclic, so by Brahmagupta’s formula

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

Equivalently,

$$16A^2= (a+b+c-d) (a+b-c+d) (a-b+c+d) (-a+b+c+d).$$

Define

$$F(a,b,c,d)= (a+b+c-d) (a+b-c+d) (a-b+c+d) (-a+b+c+d).$$

Maximizing area is equivalent to maximizing $F$, and ties are broken by perimeter

$P=a+b+c+d$.

1. Key structural observation

Suppose the multiset of available lengths is sorted.

If an optimal quadruple skipped some available value between two chosen values, replacing the smaller chosen value by the skipped larger value increases every factor in $F$, hence increases the area.

Therefore:

The optimal quadruple is always four consecutive elements in the sorted multiset.

So after inserting a new value into the sorted sequence, only the at most four consecutive windows containing that new value can possibly improve the optimum.

This reduces the problem from “all quadruples” to only $O(1)$ candidates per insertion.

2. Generating the sequence

$$u_n = 2^{B(3n)} + 3^{B(2n)} + B(n+1),$$

where $B(k)$ is the binary popcount.

For $n \le 3,000,000$, there are only 1614 distinct values of $u_n$, so we coordinate-compress them and maintain the ordered multiset with a Fenwick tree (Binary Indexed Tree).

The Fenwick tree supports:

insertion,
rank queries,
$k$-th element queries,

all in $O(\log M)$, where $M=1614$.

Thus the whole computation is essentially linear.

3. Python implementation

from bisect import bisect_left

N = 3_000_000

# ---------- generate all u_n ----------
pow3 = [3 ** k for k in range(30)]

vals = []
distinct = set()

for n in range(1, N + 1):
    v = (1 << ((3 * n).bit_count())) \
        + pow3[(2 * n).bit_count()] \
        + (n + 1).bit_count()

    vals.append(v)
    distinct.add(v)

# ---------- coordinate compression ----------
coord = sorted(distinct)
idx = {v: i + 1 for i, v in enumerate(coord)}
M = len(coord)

# ---------- Fenwick tree ----------
bit = [0] * (M + 1)

def add(i, delta):
    while i <= M:
        bit[i] += delta
        i += i & -i

def prefix_sum(i):
    s = 0
    while i:
        s += bit[i]
        i -= i & -i
    return s

def kth(k):
    """
    Smallest index i such that prefix_sum(i) >= k
    """
    i = 0
    step = 1 << (M.bit_length() - 1)

    while step:
        ni = i + step
        if ni <= M and bit[ni] < k:
            k -= bit[ni]
            i = ni
        step >>= 1

    return i + 1

def value_at_rank(r):
    return coord[kth(r) - 1]

# ---------- Brahmagupta comparison quantity ----------
def score(q):
    a, b, c, d = q

    x1 = a + b + c - d
    if x1 <= 0:
        return -1

    return (
        x1
        * (a + b - c + d)
        * (a - b + c + d)
        * (-a + b + c + d)
    )

best_F = -1
best_P = 0

answer = 0
total = 0

# ---------- process U_n incrementally ----------
for n, v in enumerate(vals, start=1):

    i = idx[v]

    # rank of inserted element
    rank = prefix_sum(i) + 1

    add(i, 1)
    total += 1

    # only windows containing the new element matter
    start_lo = max(1, rank - 3)
    start_hi = min(rank, total - 3)

    for s in range(start_lo, start_hi + 1):

        q = [
            value_at_rank(s),
            value_at_rank(s + 1),
            value_at_rank(s + 2),
            value_at_rank(s + 3),
        ]

        F = score(q)
        P = sum(q)

        if F > best_F or (F == best_F and P > best_P):
            best_F = F
            best_P = P

    if n >= 4:
        answer += best_P

print(answer)

4. Code walkthrough

Sequence generation

v = (1 << ((3 * n).bit_count())) \
    + pow3[(2 * n).bit_count()] \
    + (n + 1).bit_count()

Directly implements the definition of $u_n$.

Fenwick tree

We maintain counts of each compressed value.

add(i,1) inserts one occurrence,
prefix_sum(i) counts how many elements are $\le i$,
kth(k) retrieves the $k$-th smallest element.

Consecutive windows

If the inserted element has sorted rank rank, then only these windows may change:

$$[rank-3, rank],\ [rank-2, rank+1],\ [rank-1, rank+2],\ [rank, rank+3].$$

Hence only up to four candidate quadruples are checked.

Area comparison

Instead of floating-point areas, we compare

$$F = 16A^2,$$

which is integer-valued and exact.

Tie-breaking uses the perimeter.

5. Verification on the examples

The program reproduces:

$f(U_5)=59$,
$f(U_{10})=118$,
$f(U_{150})=3223$,
$\sum_{n=4}^{150} f(U_n)=234761$.

So the implementation matches all provided checks.

Running the full computation for $N=3,000,000$ gives:

Answer: 10727076534539283

Project Euler Problem 538