Problem 543: Prime-Sum Numbers

Project Euler Problem 543 Solution

EulerSolve provides an optimized solution for Project Euler Problem 543, Prime-Sum Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For integers $n\ge 1$ and $k\ge 1$, let $P(n,k)=1$ if $n$ can be written as a sum of exactly $k$ primes, and let $P(n,k)=0$ otherwise. Then $$S(n)=\sum_{i=1}^{n}\sum_{k=1}^{i} P(i,k).$$ Prime repetitions are allowed, but only existence matters: even if a number has many different decompositions into $k$ primes, the pair $(n,k)$ contributes only $1$. The goal is to evaluate $$\sum_{r=3}^{44} S(F_r),$$ where $F_1=1$, $F_2=1$, and $F_r=F_{r-1}+F_{r-2}$. A direct search through all prime sums is far too expensive, so the solution turns the existence questions into a closed formula for $S(n)$. Mathematical Approach The implementations split the count into the cases $k=1$, $k=2$, and $k\ge 3$. The first is pure primality, the second is a parity argument combined with the Goldbach pattern used by the solution, and the third reduces to a simple inequality. Step 1: Count admissible pairs instead of decompositions For a fixed integer $i$, the inner sum $\sum_k P(i,k)$ counts how many lengths $k$ are possible, not how many different prime partitions exist. Therefore $S(n)$ counts admissible pairs $(i,k)$ with $1\le i\le n$. This viewpoint is decisive because the implementation never enumerates actual prime sums. It only determines whether a representation of length $k$ exists....

Detailed mathematical approach

Problem Summary

For integers $n\ge 1$ and $k\ge 1$, let $P(n,k)=1$ if $n$ can be written as a sum of exactly $k$ primes, and let $P(n,k)=0$ otherwise. Then

$$S(n)=\sum_{i=1}^{n}\sum_{k=1}^{i} P(i,k).$$

Prime repetitions are allowed, but only existence matters: even if a number has many different decompositions into $k$ primes, the pair $(n,k)$ contributes only $1$. The goal is to evaluate

$$\sum_{r=3}^{44} S(F_r),$$

where $F_1=1$, $F_2=1$, and $F_r=F_{r-1}+F_{r-2}$. A direct search through all prime sums is far too expensive, so the solution turns the existence questions into a closed formula for $S(n)$.

Mathematical Approach

The implementations split the count into the cases $k=1$, $k=2$, and $k\ge 3$. The first is pure primality, the second is a parity argument combined with the Goldbach pattern used by the solution, and the third reduces to a simple inequality.

Step 1: Count admissible pairs instead of decompositions

For a fixed integer $i$, the inner sum $\sum_k P(i,k)$ counts how many lengths $k$ are possible, not how many different prime partitions exist. Therefore $S(n)$ counts admissible pairs $(i,k)$ with $1\le i\le n$.

This viewpoint is decisive because the implementation never enumerates actual prime sums. It only determines whether a representation of length $k$ exists.

Step 2: Handle the cases $k=1$ and $k=2$

When $k=1$, a number contributes exactly when it is prime, so

$$A_1(n)=\pi(n).$$

When $k=2$, parity forces a split. In the numeric range needed by the problem, the implementation counts every even integer $i\ge 4$ as a valid two-prime sum, giving

$$A_{2,\mathrm{even}}(n)=\max\left(\left\lfloor\frac{n}{2}\right\rfloor-1,0\right).$$

An odd sum of two primes must be $2+p$ with $p$ an odd prime, so the odd contribution is

$$A_{2,\mathrm{odd}}(n)=\max\left(\pi(n-2)-1,0\right).$$

Step 3: Show that for $k\ge 3$, only the bound $n\ge 2k$ matters

If $n$ is written as a sum of $k$ primes, then certainly $n\ge 2k$, because every prime is at least $2$.

Conversely, for the range handled by the solution, $n\ge 2k$ is also sufficient once $k\ge 3$. Remove $k-3$ copies of $2$, leaving the remainder

$$n-2(k-3).$$

If this remainder is odd, then it is at least $7$, so it can be expressed as a sum of three primes by the weak Goldbach theorem. If it is even, then it is at least $6$; write it as $2$ plus an even number at least $4$, and the remaining even part is handled by the same Goldbach-type two-prime case used above.

Therefore, for all $k\ge 3$,

$$P(n,k)=1 \iff n\ge 2k.$$

Step 4: Sum all contributions with at least three primes

Let $i=2t$ be even. Then the admissible lengths are $k=3,4,\dots,t$, so the number of valid $k$ values is $t-2$ whenever $t\ge 3$. Summing over all even $i\le n$ gives

$$A_{\ge 3,\mathrm{even}}(n)=\sum_{t=3}^{\lfloor n/2\rfloor}(t-2)=\max\left(\frac{(m-1)(m-2)}{2},0\right),\qquad m=\left\lfloor\frac{n}{2}\right\rfloor.$$

Now let $i=2t+1$ be odd. The admissible lengths are again $k=3,4,\dots,t$, so the count is $t-2$ for $t\ge 3$. Hence

$$A_{\ge 3,\mathrm{odd}}(n)=\sum_{t=3}^{\lfloor (n-1)/2\rfloor}(t-2)=\max\left(\frac{(m'-1)(m'-2)}{2},0\right),\qquad m'=\left\lfloor\frac{n-1}{2}\right\rfloor.$$

Step 5: Combine the pieces into one closed formula

Adding the one-prime, two-prime, and at-least-three-prime parts yields

$$\boxed{S(n)=\pi(n)+\max\left(\left\lfloor\frac{n}{2}\right\rfloor-1,0\right)+\max\left(\pi(n-2)-1,0\right)+\max\left(\frac{(m-1)(m-2)}{2},0\right)+\max\left(\frac{(m'-1)(m'-2)}{2},0\right),}$$

with

$$m=\left\lfloor\frac{n}{2}\right\rfloor,\qquad m'=\left\lfloor\frac{n-1}{2}\right\rfloor.$$

This is exactly the expression evaluated by the implementations.

Worked Example: $S(10)=20$

The sample checkpoint follows immediately from the decomposition above.

First, $\pi(10)=4$, corresponding to $2,3,5,7$.

For $k=2$, the even values are $4,6,8,10$, contributing $4$, and the odd values are $5,7,9$, contributing $3$.

For $k\ge 3$, the even numbers contribute $1+2+3=6$, because $6$ allows only $k=3$, $8$ allows $k=3,4$, and $10$ allows $k=3,4,5$.

The odd numbers contribute $1+2=3$, because $7$ allows only $k=3$ and $9$ allows $k=3,4$.

Therefore

$$S(10)=4+4+3+6+3=20.$$

Step 6: Apply the formula to the Fibonacci inputs

Once $S(n)$ is available in closed form, the outer task is simple: generate $F_3,F_4,\dots,F_{44}$ iteratively and add the corresponding values. No search over prime sums is needed any more; each term reduces to two prime-counting queries and a few integer formulas.

How the Code Works

The C++, Python, and Java implementations first build a table of primes and a prefix table for $\pi(x)$ up to five million. Larger prime-counting queries are handled by a Lehmer-style recursion with memoization, using integer square-root, cube-root, and fourth-root cutoffs to break the problem into smaller subproblems.

For each Fibonacci value $F_r$, the implementation asks only for $\pi(F_r)$ and $\pi(F_r-2)$. Those two counts are inserted into the closed formula above, and the $k=2$ and $k\ge 3$ terms are computed directly with integer arithmetic. The final loop accumulates these values for $r=3$ through $44$. The same logic also reproduces the checkpoints $S(10)=20$, $S(100)=2402$, and $S(1000)=248838$.

Complexity Analysis

Let $B=5{,}000{,}000$ be the sieve limit. Building the initial prime table costs $O(B\log\log B)$ time and $O(B)$ memory. After that, each evaluation of $S(n)$ performs two sublinear prime-counting queries plus $O(1)$ arithmetic. Since the outer sum contains only $42$ Fibonacci arguments, the total runtime is dominated by the prime-counting routine, while the closed-form part is negligible by comparison.

Footnotes and References

Problem page: https://projecteuler.net/problem=543
Prime-counting function: Wikipedia — Prime-counting function
Goldbach's conjecture: Wikipedia — Goldbach's conjecture
Weak Goldbach theorem: Wikipedia — Goldbach's weak conjecture
Fibonacci number: Wikipedia — Fibonacci number

Problem 543 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;
using u128 = unsigned __int128;

static u64 isqrt_u64(u64 n) {
    u64 x = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((x + 1ULL) <= n / (x + 1ULL)) ++x;
    while (x > 0ULL && x > n / x) --x;
    return x;
}

static u64 icbrt_u64(u64 n) {
    u64 x = static_cast<u64>(std::cbrt(static_cast<long double>(n)));
    while ((x + 1ULL) <= n / ((x + 1ULL) * (x + 1ULL))) ++x;
    while (x > 0ULL && x > n / (x * x)) --x;
    return x;
}

static u64 iroot4_u64(u64 n) {
    const long double x = static_cast<long double>(n);
    u64 r = static_cast<u64>(std::sqrt(std::sqrt(x)));
    auto pow4 = [](u64 y) -> u128 { return static_cast<u128>(y) * y * y * y; };
    while (pow4(r + 1ULL) <= n) ++r;
    while (r > 0ULL && pow4(r) > n) --r;
    return r;
}

class PrimeCounter {
public:
    explicit PrimeCounter(int sieve_limit = 5'000'000) { build_sieve(sieve_limit); }

    u64 pi(u64 n) {
        if (n <= static_cast<u64>(sieve_limit_)) {
            return pi_small_[static_cast<std::size_t>(n)];
        }
        auto it = pi_cache_.find(n);
        if (it != pi_cache_.end()) return it->second;

        const u64 a = pi(iroot4_u64(n));
        const u64 b = pi(isqrt_u64(n));
        const u64 c = pi(icbrt_u64(n));

        i64 sum = static_cast<i64>(phi(n, static_cast<int>(a))) +
                  static_cast<i64>((b + a - 2ULL) * (b - a + 1ULL) / 2ULL);

        for (u64 i = a + 1ULL; i <= b; ++i) {
            const u64 w = n / static_cast<u64>(primes_[static_cast<std::size_t>(i - 1ULL)]);
            sum -= static_cast<i64>(pi(w));
            if (i <= c) {
                const u64 lim = pi(isqrt_u64(w));
                for (u64 j = i; j <= lim; ++j) {
                    const u64 pj = static_cast<u64>(primes_[static_cast<std::size_t>(j - 1ULL)]);
                    sum -= static_cast<i64>(pi(w / pj) - (j - 1ULL));
                }
            }
        }

        const u64 out = static_cast<u64>(sum);
        pi_cache_[n] = out;
        return out;
    }

private:
    void build_sieve(int limit) {
        sieve_limit_ = limit;
        std::vector<bool> is_comp(static_cast<std::size_t>(limit + 1), false);
        pi_small_.assign(static_cast<std::size_t>(limit + 1), 0ULL);
        for (int i = 2; i <= limit; ++i) {
            if (!is_comp[static_cast<std::size_t>(i)]) {
                primes_.push_back(i);
                if (i <= limit / i) {
                    for (int j = i * i; j <= limit; j += i) is_comp[static_cast<std::size_t>(j)] = true;
                }
            }
            pi_small_[static_cast<std::size_t>(i)] =
                pi_small_[static_cast<std::size_t>(i - 1)] + (!is_comp[static_cast<std::size_t>(i)] ? 1ULL : 0ULL);
        }
    }

    u64 phi(u64 x, int s) {
        if (s == 0) return x;
        if (s == 1) return x - x / 2ULL;
        if (s == 2) return x - x / 2ULL - x / 3ULL + x / 6ULL;
        if (s == 3) return x - x / 2ULL - x / 3ULL - x / 5ULL + x / 6ULL + x / 10ULL + x / 15ULL - x / 30ULL;
        if (x <= static_cast<u64>(sieve_limit_) && static_cast<u64>(primes_[static_cast<std::size_t>(s - 1)]) >= x) return 1ULL;

        // (x,s) packing: safe here because s <= pi(iroot4(n)) <= 64 for our n-range.
        const u64 key = (x << 6U) ^ static_cast<u64>(s);
        auto it = phi_cache_.find(key);
        if (it != phi_cache_.end()) return it->second;

        const u64 ps = static_cast<u64>(primes_[static_cast<std::size_t>(s - 1)]);
        const u64 out = phi(x, s - 1) - phi(x / ps, s - 1);
        phi_cache_[key] = out;
        return out;
    }

    int sieve_limit_ = 0;
    std::vector<int> primes_;
    std::vector<u64> pi_small_;
    std::unordered_map<u64, u64> pi_cache_;
    std::unordered_map<u64, u64> phi_cache_;
};

static u128 S_of(u64 n, PrimeCounter& pc) {
    const u64 pi_n = pc.pi(n);
    const u64 pi_n_minus2 = (n >= 2 ? pc.pi(n - 2) : 0ULL);

    const u64 m = n / 2ULL;
    const u64 even_2prime = (m >= 2 ? (m - 1ULL) : 0ULL);  // even i>=4
    const u64 odd_2prime = (pi_n_minus2 >= 1 ? (pi_n_minus2 - 1ULL) : 0ULL);  // odd i=2+p, p odd prime

    const u128 sum_even_kge3 = (m >= 3 ? (static_cast<u128>(m - 1ULL) * (m - 2ULL)) / 2U : 0U);
    const u64 m2 = (n >= 1 ? (n - 1ULL) / 2ULL : 0ULL);
    const u128 sum_odd_kge3 = (m2 >= 3 ? (static_cast<u128>(m2 - 1ULL) * (m2 - 2ULL)) / 2U : 0U);

    return static_cast<u128>(pi_n) + static_cast<u128>(even_2prime) + static_cast<u128>(odd_2prime) + sum_even_kge3 +
           sum_odd_kge3;
}

static std::string to_string_u128(u128 value) {
    if (value == 0) return "0";
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    PrimeCounter pc;

    // Validation points from the problem statement.
    assert(S_of(10, pc) == 20U);
    assert(S_of(100, pc) == 2402U);
    assert(S_of(1000, pc) == 248838U);

    u64 f0 = 0, f1 = 1;
    u128 total = 0;
    for (int k = 2; k <= 44; ++k) {
        const u64 f = f0 + f1;
        f0 = f1;
        f1 = f;
        if (k >= 3) total += S_of(f, pc);
    }

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

Python

import math
import sys

sys.setrecursionlimit(20000)

def isqrt_u64(n):
    if n == 0: return 0
    return math.isqrt(n)

def icbrt_u64(n):
    if n == 0: return 0
    r = int(math.floor(n**(1/3.0)))
    while (r + 1)**3 <= n:
        r += 1
    while r**3 > n:
        r -= 1
    return r

def iroot4_u64(n):
    if n == 0: return 0
    r = int(math.floor(n**(1/4.0)))
    while (r + 1)**4 <= n:
        r += 1
    while r**4 > n:
        r -= 1
    return r

class PrimeCounter:
    def __init__(self, sieve_limit=5000000):
        self.sieve_limit = sieve_limit
        is_comp = [False] * (sieve_limit + 1)
        self.pi_small = [0] * (sieve_limit + 1)
        self.primes = []
        
        for i in range(2, sieve_limit + 1):
            if not is_comp[i]:
                self.primes.append(i)
                for j in range(i * i, sieve_limit + 1, i):
                    is_comp[j] = True
            self.pi_small[i] = self.pi_small[i - 1] + (1 if not is_comp[i] else 0)
            
        self.pi_cache = {}
        self.phi_cache = {}
        
    def pi(self, n):
        if n <= self.sieve_limit:
            return self.pi_small[n]
        if n in self.pi_cache:
            return self.pi_cache[n]
            
        a = self.pi(iroot4_u64(n))
        b = self.pi(isqrt_u64(n))
        c = self.pi(icbrt_u64(n))
        
        ans = self.phi(n, a) + (b + a - 2) * (b - a + 1) // 2
        
        for i in range(a + 1, b + 1):
            w = n // self.primes[i - 1]
            ans -= self.pi(w)
            if i <= c:
                lim = self.pi(isqrt_u64(w))
                for j in range(i, lim + 1):
                    pj = self.primes[j - 1]
                    ans -= self.pi(w // pj) - (j - 1)
                    
        self.pi_cache[n] = ans
        return ans
        
    def phi(self, x, s):
        if s == 0: return x
        if s == 1: return x - x // 2
        if s == 2: return x - x // 2 - x // 3 + x // 6
        if s == 3: return x - x // 2 - x // 3 - x // 5 + x // 6 + x // 10 + x // 15 - x // 30
        if x <= self.sieve_limit and self.primes[s - 1] >= x:
            return 1
            
        key = (x << 6) ^ s
        if key in self.phi_cache:
            return self.phi_cache[key]
            
        ps = self.primes[s - 1]
        ans = self.phi(x, s - 1) - self.phi(x // ps, s - 1)
        self.phi_cache[key] = ans
        return ans

def S_of(n, pc):
    pi_n = pc.pi(n)
    pi_n_minus2 = pc.pi(n - 2) if n >= 2 else 0
    
    m = n // 2
    even_2prime = m - 1 if m >= 2 else 0
    odd_2prime = pi_n_minus2 - 1 if pi_n_minus2 >= 1 else 0
    
    sum_even_kge3 = (m - 1) * (m - 2) // 2 if m >= 3 else 0
    m2 = (n - 1) // 2 if n >= 1 else 0
    sum_odd_kge3 = (m2 - 1) * (m2 - 2) // 2 if m2 >= 3 else 0
    
    return pi_n + even_2prime + odd_2prime + sum_even_kge3 + sum_odd_kge3

def solve():
    pc = PrimeCounter(5000000)
    f0 = 0
    f1 = 1
    total = 0
    for k in range(2, 45):
        f = f0 + f1
        f0 = f1
        f1 = f
        if k >= 3:
            total += S_of(f, pc)
    return str(total)

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

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler543 {
    static long isqrt(long n) {
        if (n == 0)
            return 0;
        long x = (long) Math.sqrt((double) n);
        while ((x + 1) <= n / (x + 1))
            x++;
        while (x > 0 && x > n / x)
            x--;
        return x;
    }

    static long icbrt(long n) {
        if (n == 0)
            return 0;
        long x = (long) Math.cbrt((double) n);
        while ((x + 1) <= n / ((x + 1) * (x + 1)))
            x++;
        while (x > 0 && x > n / (x * x))
            x--;
        return x;
    }

    static long iroot4(long n) {
        if (n == 0)
            return 0;
        long r = (long) Math.sqrt(Math.sqrt((double) n));
        while (Math.pow(r + 1, 4) <= n)
            r++;
        while (r > 0 && Math.pow(r, 4) > n)
            r--;
        return r;
    }

    static class PrimeCounter {
        int sieveLimit;
        List<Integer> primes;
        long[] piSmall;
        Map<Long, Long> piCache;
        Map<Long, Long> phiCache;

        PrimeCounter(int sieveLimit) {
            this.sieveLimit = sieveLimit;
            boolean[] isComp = new boolean[sieveLimit + 1];
            piSmall = new long[sieveLimit + 1];
            primes = new ArrayList<>();

            for (int i = 2; i <= sieveLimit; i++) {
                if (!isComp[i]) {
                    primes.add(i);
                    if (i <= sieveLimit / i) {
                        for (int j = i * i; j <= sieveLimit; j += i) {
                            isComp[j] = true;
                        }
                    }
                }
                piSmall[i] = piSmall[i - 1] + (!isComp[i] ? 1 : 0);
            }
            piCache = new HashMap<>();
            phiCache = new HashMap<>();
        }

        long pi(long n) {
            if (n <= sieveLimit) {
                return piSmall[(int) n];
            }
            Long cached = piCache.get(n);
            if (cached != null)
                return cached;

            long a = pi(iroot4(n));
            long b = pi(isqrt(n));
            long c = pi(icbrt(n));

            long sum = phi(n, (int) a) + (b + a - 2) * (b - a + 1) / 2;

            for (long i = a + 1; i <= b; i++) {
                long w = n / primes.get((int) (i - 1));
                sum -= pi(w);
                if (i <= c) {
                    long lim = pi(isqrt(w));
                    for (long j = i; j <= lim; j++) {
                        long pj = primes.get((int) (j - 1));
                        sum -= pi(w / pj) - (j - 1);
                    }
                }
            }

            piCache.put(n, sum);
            return sum;
        }

        long phi(long x, int s) {
            if (s == 0)
                return x;
            if (s == 1)
                return x - x / 2;
            if (s == 2)
                return x - x / 2 - x / 3 + x / 6;
            if (s == 3)
                return x - x / 2 - x / 3 - x / 5 + x / 6 + x / 10 + x / 15 - x / 30;
            if (x <= sieveLimit && primes.get(s - 1) >= x)
                return 1;

            long key = (x << 6) ^ s;
            Long cached = phiCache.get(key);
            if (cached != null)
                return cached;

            long ps = primes.get(s - 1);
            long ans = phi(x, s - 1) - phi(x / ps, s - 1);
            phiCache.put(key, ans);
            return ans;
        }
    }

    static long S_of(long n, PrimeCounter pc) {
        long pi_n = pc.pi(n);
        long pi_n_minus2 = n >= 2 ? pc.pi(n - 2) : 0;

        long m = n / 2;
        long even_2prime = m >= 2 ? (m - 1) : 0;
        long odd_2prime = pi_n_minus2 >= 1 ? (pi_n_minus2 - 1) : 0;

        long sum_even_kge3 = m >= 3 ? (m - 1) * (m - 2) / 2 : 0;
        long m2 = n >= 1 ? (n - 1) / 2 : 0;
        long sum_odd_kge3 = m2 >= 3 ? (m2 - 1) * (m2 - 2) / 2 : 0;

        return pi_n + even_2prime + odd_2prime + sum_even_kge3 + sum_odd_kge3;
    }

    public static String solve() {
        PrimeCounter pc = new PrimeCounter(5000000);
        long f0 = 0;
        long f1 = 1;
        long total = 0;
        for (int k = 2; k <= 44; k++) {
            long f = f0 + f1;
            f0 = f1;
            f1 = f;
            if (k >= 3) {
                total += S_of(f, pc);
            }
        }
        return Long.toString(total);
    }

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