Problem 753: Fermat Equation

View on Project Euler

Project Euler Problem 753 Solution

EulerSolve provides an optimized solution for Project Euler Problem 753, Fermat Equation, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a prime \(p\), define $$F(p)=\#\left\{(a,b,c)\in\{1,\dots,p-1\}^3 : a^3+b^3\equiv c^3 \pmod p\right\}.$$ So \(F(p)\) counts ordered triples of nonzero residues modulo \(p\) that satisfy the Fermat-type congruence \(a^3+b^3=c^3\). The goal is to compute $$\sum_{\substack{p\lt 6\cdot 10^6 \\ p\text{ prime}}} F(p).$$ A direct \(O(p^3)\) count for every prime would be far too slow, so the solution uses the structure of the cube map on \(\mathbb F_p^\times\) and a closed formula for the difficult congruence class \(p\equiv 1\pmod 3\). Mathematical Approach The key observation is that the number of cube roots of a nonzero residue depends entirely on the residue class of \(p\) modulo \(3\). That turns the whole problem into a clean case split. Step 1: Rewrite the Counting Problem in \(\mathbb F_p^\times\) Because \(a,b,c\) are chosen from \(\{1,\dots,p-1\}\), all three variables live in the multiplicative group \(\mathbb F_p^\times\), which has size \(p-1\). For a nonzero residue \(s\), let \(N_p(s)\) be the number of nonzero solutions of \(x^3=s\). Then $$F(p)=\sum_{a\in\mathbb F_p^\times}\sum_{b\in\mathbb F_p^\times} N_p(a^3+b^3).$$ So the job is to understand how many cube roots a residue has and how often \(a^3+b^3\) lands in each class. Step 2: Handle the Easy Branch \(p\equiv 2\pmod 3\) The group \(\mathbb F_p^\times\) is cyclic of size \(p-1\)....

Detailed mathematical approach

Problem Summary

For a prime \(p\), define

$$F(p)=\#\left\{(a,b,c)\in\{1,\dots,p-1\}^3 : a^3+b^3\equiv c^3 \pmod p\right\}.$$

So \(F(p)\) counts ordered triples of nonzero residues modulo \(p\) that satisfy the Fermat-type congruence \(a^3+b^3=c^3\). The goal is to compute

$$\sum_{\substack{p\lt 6\cdot 10^6 \\ p\text{ prime}}} F(p).$$

A direct \(O(p^3)\) count for every prime would be far too slow, so the solution uses the structure of the cube map on \(\mathbb F_p^\times\) and a closed formula for the difficult congruence class \(p\equiv 1\pmod 3\).

Mathematical Approach

The key observation is that the number of cube roots of a nonzero residue depends entirely on the residue class of \(p\) modulo \(3\). That turns the whole problem into a clean case split.

Step 1: Rewrite the Counting Problem in \(\mathbb F_p^\times\)

Because \(a,b,c\) are chosen from \(\{1,\dots,p-1\}\), all three variables live in the multiplicative group \(\mathbb F_p^\times\), which has size \(p-1\).

For a nonzero residue \(s\), let \(N_p(s)\) be the number of nonzero solutions of \(x^3=s\). Then

$$F(p)=\sum_{a\in\mathbb F_p^\times}\sum_{b\in\mathbb F_p^\times} N_p(a^3+b^3).$$

So the job is to understand how many cube roots a residue has and how often \(a^3+b^3\) lands in each class.

Step 2: Handle the Easy Branch \(p\equiv 2\pmod 3\)

The group \(\mathbb F_p^\times\) is cyclic of size \(p-1\). If \(p\equiv 2\pmod 3\), then \(\gcd(3,p-1)=1\), so the map

$$x\longmapsto x^3$$

is a permutation of \(\mathbb F_p^\times\). Therefore every nonzero residue has exactly one nonzero cube root.

There are \((p-1)^2\) ordered pairs \((a,b)\). Among them, exactly \(p-1\) satisfy \(a^3+b^3\equiv 0\pmod p\), because for each nonzero \(a\) there is exactly one choice of \(b\), namely \(b\equiv -a\), that cancels the sum.

Those pairs contribute nothing, since \(c\) is required to be nonzero. Every other pair contributes exactly one value of \(c\). Hence

$$F(p)=(p-1)^2-(p-1)=(p-1)(p-2)\qquad (p\equiv 2\pmod 3).$$

The small primes fit naturally into this picture: \(F(2)=0\), while \(F(3)=2\) can be checked directly.

Step 3: Understand the Hard Branch \(p\equiv 1\pmod 3\)

If \(p\equiv 1\pmod 3\), then \(3\mid (p-1)\). The cube map is no longer bijective: its kernel has size \(3\), so every nonzero cubic residue has exactly three cube roots, and non-cubic residues have none.

What remains is to measure how often \(a^3+b^3\) lands in the cubic-residue set. That is the nontrivial part of the problem, and the implementations use the classical closed form

$$F(p)=(p-1)(p-8+t)\qquad (p\equiv 1\pmod 3),$$

where the integer \(t\) comes from a quadratic representation of \(p\). This formula is the cubic-character/Jacobi-sum evaluation encoded by the program.

Step 4: Determine \(t\) from the Representation \(4p=27u^2+t^2\)

For every prime \(p\equiv 1\pmod 3\), there exist integers \(u\) and \(t\) such that

$$4p=27u^2+t^2,\qquad t\equiv 1\pmod 3.$$

Up to signs, this representation is unique. The sign of \(t\) matters in the formula, so once \(|t|\) is found the correct sign is chosen by the congruence condition \(t\equiv 1\pmod 3\).

Substituting that value into the closed form immediately gives \(F(p)\).

Step 5: Assemble the Prime Sum

The complete evaluation is therefore

$$F(p)= \begin{cases} 0,&p=2,\\ 2,&p=3,\\ (p-1)(p-2),&p\equiv 2\pmod 3,\\ (p-1)(p-8+t),&p\equiv 1\pmod 3,\ 4p=27u^2+t^2,\ t\equiv 1\pmod 3. \end{cases}$$

After applying this formula prime by prime, the required answer is obtained by summing over all primes below \(6\cdot 10^6\).

Worked Example

Two small primes show both branches.

For \(p=5\), we are in the easy branch \(p\equiv 2\pmod 3\), so

$$F(5)=(5-1)(5-2)=4\cdot 3=12.$$

For \(p=19\), we are in the hard branch \(p\equiv 1\pmod 3\). We have

$$4p=76=27\cdot 1^2+7^2,$$

and \(7\equiv 1\pmod 3\), so \(t=7\). Therefore

$$F(19)=(19-1)(19-8+7)=18\cdot 18=324.$$

These values agree with the formula used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same plan. First they build a prime sieve up to \(6\cdot 10^6\), so each prime can be processed exactly once.

Next they precompute the monotone list

$$27\cdot 1^2,\ 27\cdot 2^2,\ 27\cdot 3^2,\ \dots$$

up to \(4\cdot 6\cdot 10^6\). This turns the search for the representation \(4p=27u^2+t^2\) into a scan over candidate values of \(27u^2\).

For each prime \(p\):

If \(p=3\), the contribution is \(2\).

If \(p\equiv 2\pmod 3\), the implementation adds \((p-1)(p-2)\).

If \(p\equiv 1\pmod 3\), it loops through the precomputed \(27u^2\) values not exceeding \(4p\), computes

$$4p-27u^2,$$

tests whether that remainder is a perfect square, and when a square is found uses its root as \(|t|\). The sign is then fixed so that \(t\equiv 1\pmod 3\), and the closed formula \((p-1)(p-8+t)\) is added to the total.

All three implementations accumulate the final sum in 64-bit arithmetic. One implementation also checks the formula against direct enumeration for small primes before printing the final total.

Complexity Analysis

Let \(L=6\cdot 10^6\). The prime sieve costs \(O(L\log\log L)\) time and \(O(L)\) memory. Precomputing the values \(27u^2\) takes \(O(\sqrt L)\) time and memory.

For primes with \(p\equiv 1\pmod 3\), the representation search scans \(u\) while \(27u^2\le 4p\), so a conservative upper bound is \(O(\sqrt p)\) work per such prime. Summed over all primes, this gives the practical estimate

$$O\!\left(L\log\log L+\pi(L)\sqrt L\right)$$

with \(O(L)\) total memory. In practice it is faster than this bound suggests, because the search stops as soon as the first valid representation is found.

Footnotes and References

  1. Problem page: Project Euler 753
  2. Finite fields: Wikipedia — Finite field
  3. Cubic residues: Wikipedia — Cubic residue
  4. Jacobi sums: Wikipedia — Jacobi sum
  5. Binary quadratic forms: Wikipedia — Binary quadratic form

Problem 753 source code

C++

#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using i64 = std::int64_t;

i64 isqrt_i64(const i64 n) {
    long double d = static_cast<long double>(n);
    i64 x = static_cast<i64>(std::sqrt(d));
    while ((x + 1) * (x + 1) <= n) {
        ++x;
    }
    while (x * x > n) {
        --x;
    }
    return x;
}

i64 F_formula(const int p, const std::vector<i64>& twenty_seven_u2) {
    if (p == 3) {
        return 2;
    }
    if (p < 5 || p % 3 == 0) {
        return 0;
    }
    if (p % 3 == 2) {
        return static_cast<i64>(p - 1) * static_cast<i64>(p - 2);
    }

    const i64 target = 4LL * p;
    i64 t_sel = 0;
    for (const i64 v : twenty_seven_u2) {
        if (v > target) {
            break;
        }
        const i64 rem = target - v;
        const i64 t = isqrt_i64(rem);
        if (t * t == rem) {
            t_sel = (t % 3 == 1) ? t : -t;
            break;
        }
    }
    assert(t_sel != 0);

    return static_cast<i64>(p - 1) * static_cast<i64>(p - 8 + t_sel);
}

i64 F_bruteforce(const int p) {
    i64 count = 0;
    std::vector<int> cubes(static_cast<std::size_t>(p), 0);
    for (int x = 0; x < p; ++x) {
        cubes[static_cast<std::size_t>(x)] = static_cast<int>((1LL * x * x * x) % p);
    }

    for (int a = 1; a < p; ++a) {
        const int a3 = cubes[static_cast<std::size_t>(a)];
        for (int b = 1; b < p; ++b) {
            const int s = (a3 + cubes[static_cast<std::size_t>(b)]) % p;
            for (int c = 1; c < p; ++c) {
                if (cubes[static_cast<std::size_t>(c)] == s) {
                    ++count;
                }
            }
        }
    }
    return count;
}

i64 solve() {
    constexpr int kLimit = 6'000'000;
    std::vector<bool> is_prime(static_cast<std::size_t>(kLimit), true);
    is_prime[0] = false;
    is_prime[1] = false;
    for (int i = 2; 1LL * i * i < kLimit; ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j < kLimit; j += i) {
            is_prime[static_cast<std::size_t>(j)] = false;
        }
    }

    std::vector<i64> twenty_seven_u2;
    for (i64 u = 1;; ++u) {
        const i64 v = 27LL * u * u;
        if (v > 4LL * kLimit) {
            break;
        }
        twenty_seven_u2.push_back(v);
    }

    i64 total = 0;
    for (int p = 2; p < kLimit; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        total += F_formula(p, twenty_seven_u2);
    }

    return total;
}

}  // namespace

int main() {
    std::vector<i64> test_u2;
    for (i64 u = 1; 27LL * u * u <= 4LL * 10'000; ++u) {
        test_u2.push_back(27LL * u * u);
    }

    assert(F_formula(5, test_u2) == 12);
    assert(F_formula(7, test_u2) == 0);

    for (int p = 5; p < 200; ++p) {
        bool prime = true;
        for (int d = 2; d * d <= p; ++d) {
            if (p % d == 0) {
                prime = false;
                break;
            }
        }
        if (!prime) {
            continue;
        }
        assert(F_formula(p, test_u2) == F_bruteforce(p));
    }

    std::cout << solve() << '\n';
    return 0;
}

Python

import math

def solve():
    LIMIT = 6000000
    is_prime = bytearray(b'\x01' * LIMIT)
    is_prime[0] = is_prime[1] = 0
    for i in range(2, int(LIMIT**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, LIMIT, i): is_prime[j] = 0

    def isqrt(n):
        r = int(math.isqrt(n))
        while (r+1)*(r+1) <= n: r += 1
        while r*r > n: r -= 1
        return r

    twu = []
    u = 1
    while 27*u*u <= 4*LIMIT:
        twu.append(27*u*u); u += 1

    total = 0
    for p in range(2, LIMIT):
        if not is_prime[p]: continue
        if p == 3: total += 2; continue
        if p < 5 or p % 3 == 0: continue
        if p % 3 == 2:
            total += (p-1)*(p-2); continue
        target = 4*p; t_sel = 0
        for v in twu:
            if v > target: break
            rem = target - v
            t = isqrt(rem)
            if t*t == rem:
                t_sel = t if t%3 == 1 else -t; break
        total += (p-1)*(p - 8 + t_sel)

    return str(total)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;

public class Euler753 {
    static long isqrtI64(long n) {
        if (n == 0)
            return 0;
        long x = (long) Math.sqrt(n);
        while ((x + 1) * (x + 1) <= n)
            x++;
        while (x * x > n)
            x--;
        return x;
    }

    static long FFormula(int p, List<Long> twentySevenU2) {
        if (p == 3) {
            return 2;
        }
        if (p < 5 || p % 3 == 0) {
            return 0;
        }
        if (p % 3 == 2) {
            return (long) (p - 1) * (p - 2);
        }

        long target = 4L * p;
        long tSel = 0;

        for (long v : twentySevenU2) {
            if (v > target) {
                break;
            }
            long rem = target - v;
            long t = isqrtI64(rem);
            if (t * t == rem) {
                tSel = (t % 3 == 1) ? t : -t;
                break;
            }
        }

        return (long) (p - 1) * (p - 8 + tSel);
    }

    public static String solve() {
        int kLimit = 6000000;
        byte[] isPrime = new byte[kLimit];
        for (int i = 2; i < kLimit; i++) {
            isPrime[i] = 1;
        }

        for (int i = 2; (long) i * i < kLimit; ++i) {
            if (isPrime[i] != 1)
                continue;
            for (int j = i * i; j < kLimit; j += i) {
                isPrime[j] = 0;
            }
        }

        List<Long> twentySevenU2 = new ArrayList<>();
        for (long u = 1;; ++u) {
            long v = 27L * u * u;
            if (v > 4L * kLimit) {
                break;
            }
            twentySevenU2.add(v);
        }

        long total = 0;
        for (int p = 2; p < kLimit; ++p) {
            if (isPrime[p] != 1)
                continue;
            total += FFormula(p, twentySevenU2);
        }

        return Long.toString(total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}