Solution

Answer: 2032447591196869022

Let

$$a^2+b^2=N,\qquad 0\le a\le b,$$

and define $S(N)$ as the sum of all admissible values of $a$.

We must compute

$$\sum S(N),$$

where $N$ runs over all squarefree integers composed only of primes $p\equiv 1\pmod 4$ with $p<150$.

The relevant primes are

$$5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149.$$

1. Mathematical analysis

1.1 Sum of two squares theorem

A prime $p\equiv 1\pmod 4$ has a unique representation

$$p=a^2+b^2$$

up to signs and order.

For example:

$$5=1^2+2^2,\qquad 13=2^2+3^2,\qquad 17=1^2+4^2.$$

Since every allowed $N$ is squarefree and only contains such primes, every admissible $N$ is itself representable as a sum of two squares.

1.2 Brahmagupta–Fibonacci identity

The key identity is

$$(a^2+b^2)(c^2+d^2) = (ac\pm bd)^2 + (ad\mp bc)^2.$$

Equivalently, using Gaussian integers,

$$(a+bi)(c+di) = (ac-bd)+(ad+bc)i.$$

Thus every time we multiply by another prime representation, each existing representation branches into two new ones.

1.3 Counting representations

Suppose

$$N=p_1p_2\cdots p_k$$

with distinct primes $p_i\equiv1\pmod4$.

Each prime contributes one Gaussian factor

$$a_i+b_i i.$$

At every multiplication we get two independent choices, so there are

$$2^{k-1}$$

distinct solutions with $0\le a\le b$.

Since there are $16$ eligible primes, brute force over all $2^{16}-1=65535$ squarefree $N$ is completely feasible.

1.4 Constructing all representations

For each prime $p$, let its canonical representation be

$$p=a^2+b^2,\qquad a\le b.$$

Starting from $(1,0)$, repeatedly combine representations using

$$(x,y)\otimes(a,b) = \left( |xa-yb|,; xb+ya \right)$$

and

$$\left( |xb-ya|,; xa+yb \right).$$

After each step we reorder coordinates so that

$$0\le a\le b.$$

The set of all such pairs gives all representations of $N$, and

$$S(N)=\sum a.$$

Finally we sum over all nonempty subsets of the 16 primes.

2. Python implementation

from math import isqrt

# All primes p < 150 with p % 4 == 1
primes = [
    5, 13, 17, 29, 37, 41, 53, 61,
    73, 89, 97, 101, 109, 113, 137, 149
]

# ------------------------------------------------------------
# Find the unique representation p = a^2 + b^2 with a <= b
# ------------------------------------------------------------

prime_reps = {}

for p in primes:
    for a in range(1, isqrt(p) + 1):
        b2 = p - a * a
        b = isqrt(b2)

        if b * b == b2 and a <= b:
            prime_reps[p] = (a, b)
            break

# ------------------------------------------------------------
# Normalize pair so that 0 <= a <= b
# ------------------------------------------------------------

def canonical(a, b):
    a = abs(a)
    b = abs(b)

    if a > b:
        a, b = b, a

    return (a, b)

# ------------------------------------------------------------
# Combine two representations using
# (a^2+b^2)(c^2+d^2)
# ------------------------------------------------------------

def combine(rep1, rep2):
    a, b = rep1
    c, d = rep2

    # Two independent combinations
    x1 = canonical(a * c - b * d, a * d + b * c)
    x2 = canonical(a * d - b * c, a * c + b * d)

    return [x1, x2]

# ------------------------------------------------------------
# Enumerate all squarefree products
# ------------------------------------------------------------

total = 0
n = len(primes)

for mask in range(1, 1 << n):

    # Current set of representations
    reps = {(1, 0)}

    for i in range(n):
        if (mask >> i) & 1:

            pair = prime_reps[primes[i]]

            new_reps = set()

            for r in reps:
                for nr in combine(r, pair):
                    new_reps.add(nr)

            reps = new_reps

    # S(N) = sum of the smaller coordinates
    total += sum(a for a, b in reps)

print(total)

3. Code walkthrough

Step 1 — Prime representations

For every prime $p\equiv1\pmod4$, we search for integers satisfying

$$p=a^2+b^2.$$

Example:

13 = 2^2 + 3^2

so we store:

prime_reps[13] = (2, 3)

Step 2 — Canonical ordering

Because

$$(a,b),\ (b,a),\ (-a,b)$$

all represent the same solution geometrically, we normalize every pair to

$$0\le a\le b.$$

This prevents duplicates.

Step 3 — Combining representations

Using

$$(a^2+b^2)(c^2+d^2) = (ac-bd)^2+(ad+bc)^2,$$

and the second sign choice,

$$(ad-bc)^2+(ac+bd)^2,$$

we generate all representations of products.

Step 4 — Enumerating all $N$

Every squarefree $N$ corresponds to a subset of the 16 primes.

We loop through all bitmasks:

for mask in range(1, 1 << 16):

and recursively build all representations.

Step 5 — Example check: $N=65$

$$65=5\cdot13.$$

The algorithm generates:

$$1^2+8^2,\qquad 4^2+7^2.$$

Therefore

$$S(65)=1+4=5,$$

matching the statement.

4. Final verification

There are only $2^{16}-1=65535$ values of $N$, so exhaustive enumeration is safe.
Every representation is generated exactly once after normalization.
The sample $S(65)=5$ matches.
The computed total is a large positive integer, consistent with the exponential growth in representation counts.

Therefore the result is:

Answer: 2032447591196869022

Project Euler Problem 273