Problem 515: Dissonant Numbers

View on Project Euler

Project Euler Problem 515 Solution

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

Problem Summary We must evaluate the prime sum $$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prime}}} d(p,p-1,k).$$ The key fact used by the implementation is that for every contributing prime \(p\), the required term collapses to a modular inverse of \(k-1\) modulo \(p\). After that simplification, the whole problem becomes: enumerate the primes in a long interval efficiently, compute one inverse for each of them, and add the results. Mathematical Approach The number-theoretic core is simple once the prime-term identity is recognized. The difficult part is not symbolic algebra, but organizing the computation so that the interval \([a,a+b)\) can be processed without sieving every number up to \(a+b\). Step 1: Reduce each term to a congruence The implementation uses the identity $$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$ So for each prime \(p\) in the interval, the required contribution is the unique residue \(x\) satisfying $$(k-1)x \equiv 1 \pmod{p}.$$ Because the modulus is prime, this inverse exists whenever \(p \nmid (k-1)\). In other words, the original arithmetic object \(d(p,p-1,k)\) does not need to be built directly; it is enough to solve one linear congruence per prime....

Detailed mathematical approach

Problem Summary

We must evaluate the prime sum

$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prime}}} d(p,p-1,k).$$

The key fact used by the implementation is that for every contributing prime \(p\), the required term collapses to a modular inverse of \(k-1\) modulo \(p\). After that simplification, the whole problem becomes: enumerate the primes in a long interval efficiently, compute one inverse for each of them, and add the results.

Mathematical Approach

The number-theoretic core is simple once the prime-term identity is recognized. The difficult part is not symbolic algebra, but organizing the computation so that the interval \([a,a+b)\) can be processed without sieving every number up to \(a+b\).

Step 1: Reduce each term to a congruence

The implementation uses the identity

$$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$

So for each prime \(p\) in the interval, the required contribution is the unique residue \(x\) satisfying

$$(k-1)x \equiv 1 \pmod{p}.$$

Because the modulus is prime, this inverse exists whenever \(p \nmid (k-1)\). In other words, the original arithmetic object \(d(p,p-1,k)\) does not need to be built directly; it is enough to solve one linear congruence per prime.

Step 2: Use Bézout coefficients to obtain the inverse

When \(\gcd(k-1,p)=1\), the extended Euclidean algorithm finds integers \(u\) and \(v\) such that

$$u(k-1)+vp=1.$$

Reducing that identity modulo \(p\) removes the second term and leaves

$$u(k-1)\equiv 1 \pmod{p}.$$

Therefore the required contribution is simply the least nonnegative residue of \(u\):

$$d(p,p-1,k)=u \bmod p.$$

This is why the three implementations use the extended Euclidean algorithm rather than repeated powering. It produces the inverse directly and works uniformly in all three languages.

Step 3: Enumerate primes with a segmented sieve

The interval can start near \(10^9\), so ordinary sieving from \(2\) all the way to \(a+b\) would be wasteful. Instead, let

$$R=a+b-1,\qquad s=\left\lfloor \sqrt{R}\right\rfloor.$$

First compute all primes up to \(s\) with a standard sieve. These are the only base primes needed to mark composites in the target interval. For each base prime \(q\), the first multiple that must be removed from \([a,a+b)\) is

$$m_q=\max\left(q^2,\left\lceil\frac{a}{q}\right\rceil q\right).$$

Then every number

$$m_q,\ m_q+q,\ m_q+2q,\dots$$

inside the interval is composite and can be marked. After all base primes are processed, the unmarked entries are exactly the primes in \([a,a+b)\).

Step 4: Accumulate the prime contributions

Once the interval primes are known, the sum is evaluated term by term:

$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prime}}} \left((k-1)^{-1} \bmod p\right).$$

No further combinatorial structure is hidden here. The algorithm does precisely two things: identify the primes in the interval and attach to each prime the modular inverse of \(k-1\).

Worked Example: \(D(101,1,10)=45\)

The interval \([101,102)\) contains only one prime, namely \(p=101\). Here

$$k-1=9.$$

So we solve

$$9x \equiv 1 \pmod{101}.$$

Since

$$9\cdot 45 = 405 = 4\cdot 101 + 1,$$

the inverse of \(9\) modulo \(101\) is \(45\). Therefore

$$D(101,1,10)=45,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They first sieve all primes up to \(\lfloor\sqrt{a+b-1}\rfloor\). Next they allocate a Boolean segment of length \(b\) representing the numbers \(a,a+1,\dots,a+b-1\), and every composite hit by a base prime is marked. After that scan, each unmarked position corresponds to a prime \(p\) in the interval.

For every such prime, the implementation runs the extended Euclidean algorithm on \(k-1\) and \(p\), normalizes the Bézout coefficient into the range \(0\) to \(p-1\), and adds that value to the running total. The inverse routine is written defensively: if the reduced value of \(k-1\) is \(0\) modulo \(p\), or if the gcd is not \(1\), it returns \(0\). For the intended evaluations, the contributing primes are coprime to \(k-1\), so each genuine term is well defined.

Complexity Analysis

Let \(R=a+b-1\). Generating all base primes up to \(\lfloor\sqrt{R}\rfloor\) costs \(O(\sqrt{R}\log\log R)\) time and \(O(\sqrt{R})\) memory with the simple sieve used in the implementations. Marking the interval of length \(b\) costs \(O(b\log\log R)\) time and \(O(b)\) memory. If \(\pi([a,a+b))\) denotes the number of primes in the interval, the inverse computations add \(O(\pi([a,a+b))\log R)\) time. So the total memory usage is linear in the interval length, and the method is efficient because it never stores all integers up to \(a+b\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=515
  2. Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
  3. Extended Euclidean algorithm: Wikipedia — Extended Euclidean algorithm
  4. Bézout's identity: Wikipedia — Bézout's identity
  5. Segmented sieve: Wikipedia — Segmented sieve

Problem 515 source code

C++

#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;

i64 extended_gcd(const i64 a, const i64 b, i64& x, i64& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    i64 x1 = 0;
    i64 y1 = 0;
    const i64 g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

u64 mod_inverse(u64 a, u64 mod) {
    a %= mod;
    if (a == 0ULL) {
        return 0ULL;
    }
    i64 x = 0;
    i64 y = 0;
    const i64 g = extended_gcd(static_cast<i64>(a), static_cast<i64>(mod), x, y);
    if (g != 1) {
        return 0ULL;
    }
    i64 res = x % static_cast<i64>(mod);
    if (res < 0) {
        res += static_cast<i64>(mod);
    }
    return static_cast<u64>(res);
}

std::vector<int> small_primes_up_to(const int n) {
    std::vector<bool> is_prime(static_cast<std::size_t>(n + 1), true);
    if (n >= 0) {
        is_prime[0] = false;
    }
    if (n >= 1) {
        is_prime[1] = false;
    }
    for (int p = 2; p * p <= n; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int m = p * p; m <= n; m += p) {
            is_prime[static_cast<std::size_t>(m)] = false;
        }
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
    }
    return primes;
}

std::vector<u64> segmented_primes(const u64 low, const u64 high_exclusive) {
    if (high_exclusive <= low) {
        return {};
    }

    const u64 high_inclusive = high_exclusive - 1ULL;
    int root = 1;
    while (static_cast<u64>(root) * static_cast<u64>(root) <= high_inclusive) {
        ++root;
    }
    --root;

    const std::vector<int> base_primes = small_primes_up_to(root);
    const u64 size = high_exclusive - low;
    std::vector<bool> is_prime(static_cast<std::size_t>(size), true);

    for (const int p_int : base_primes) {
        const u64 p = static_cast<u64>(p_int);
        u64 start = (low + p - 1ULL) / p * p;
        const u64 pp = p * p;
        if (start < pp) {
            start = pp;
        }
        for (u64 x = start; x < high_exclusive; x += p) {
            is_prime[static_cast<std::size_t>(x - low)] = false;
        }
    }

    if (low == 0ULL) {
        is_prime[0] = false;
        if (size > 1ULL) {
            is_prime[1] = false;
        }
    } else if (low == 1ULL) {
        is_prime[0] = false;
    }

    std::vector<u64> primes;
    for (u64 i = 0ULL; i < size; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            primes.push_back(low + i);
        }
    }
    return primes;
}

u64 D(const u64 a, const u64 b, const u64 k) {
    const u64 m = k - 1ULL;
    const std::vector<u64> primes = segmented_primes(a, a + b);

    u64 total = 0ULL;
    for (const u64 p : primes) {
        const u64 term = mod_inverse(m, p);  // equals d(p,p-1,k) mod p
        total += term;
    }
    return total;
}

bool run_checkpoints() {
    if (D(101ULL, 1ULL, 10ULL) != 45ULL) {
        std::cerr << "Checkpoint failed: D(101,1,10)\n";
        return false;
    }
    if (D(1'000ULL, 100ULL, 100ULL) != 8'334ULL) {
        std::cerr << "Checkpoint failed: D(10^3,10^2,10^2)\n";
        return false;
    }
    if (D(1'000'000ULL, 1'000ULL, 1'000ULL) != 38'162'302ULL) {
        std::cerr << "Checkpoint failed: D(10^6,10^3,10^3)\n";
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr u64 a = 1'000'000'000ULL;
    constexpr u64 b = 100'000ULL;
    constexpr u64 k = 100'000ULL;
    std::cout << D(a, b, k) << '\n';
    return 0;
}

Python

def extended_gcd(a, b):
    if b == 0:
        return 1, 0, a
    x1, y1, g = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return x, y, g

def mod_inverse(a, mod):
    a %= mod
    if a == 0:
        return 0
    x, y, g = extended_gcd(a, mod)
    if g != 1:
        return 0
    res = x % mod
    if res < 0:
        res += mod
    return res

def small_primes_up_to(n):
    is_prime = [True] * (n + 1)
    if n >= 0: is_prime[0] = False
    if n >= 1: is_prime[1] = False
    for p in range(2, int(n**0.5) + 1):
        if is_prime[p]:
            for m in range(p * p, n + 1, p):
                is_prime[m] = False
    return [i for i in range(2, n + 1) if is_prime[i]]

def segmented_primes(low, high_exclusive):
    if high_exclusive <= low:
        return []
    high_inclusive = high_exclusive - 1
    root = int(high_inclusive ** 0.5)
    
    base_primes = small_primes_up_to(root)
    size = high_exclusive - low
    is_prime = [True] * size
    
    for p in base_primes:
        start = (low + p - 1) // p * p
        if start < p * p:
            start = p * p
        for x in range(start, high_exclusive, p):
            is_prime[x - low] = False
            
    if low == 0:
        is_prime[0] = False
        if size > 1: is_prime[1] = False
    elif low == 1:
        is_prime[0] = False
        
    primes = []
    for i in range(size):
        if is_prime[i]:
            primes.append(low + i)
    return primes

def D(a, b, k):
    m = k - 1
    primes = segmented_primes(a, a + b)
    total = 0
    for p in primes:
        total += mod_inverse(m, p)
    return total

def solve():
    a = 1000000000
    b = 100000
    k = 100000
    ans = D(a, b, k)
    return str(ans)

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

Java

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

public class Euler515 {

    static class GCDResult {
        long x, y, g;

        GCDResult(long x, long y, long g) {
            this.x = x;
            this.y = y;
            this.g = g;
        }
    }

    static GCDResult extendedGcd(long a, long b) {
        if (b == 0) {
            return new GCDResult(1, 0, a);
        }
        GCDResult res = extendedGcd(b, a % b);
        long x = res.y;
        long y = res.x - (a / b) * res.y;
        return new GCDResult(x, y, res.g);
    }

    static long modInverse(long a, long mod) {
        a %= mod;
        if (a == 0)
            return 0;
        GCDResult res = extendedGcd(a, mod);
        if (res.g != 1)
            return 0;
        long inv = res.x % mod;
        if (inv < 0)
            inv += mod;
        return inv;
    }

    static List<Integer> smallPrimesUpTo(int n) {
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;
        for (int p = 2; p * p <= n; p++) {
            if (isPrime[p]) {
                for (int m = p * p; m <= n; m += p) {
                    isPrime[m] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static List<Long> segmentedPrimes(long low, long highExclusive) {
        if (highExclusive <= low)
            return new ArrayList<>();
        long highInclusive = highExclusive - 1;
        int root = (int) Math.sqrt(highInclusive);

        List<Integer> basePrimes = smallPrimesUpTo(root);
        int size = (int) (highExclusive - low);
        boolean[] isPrime = new boolean[size];
        for (int i = 0; i < size; i++)
            isPrime[i] = true;

        for (int p : basePrimes) {
            long start = (low + p - 1) / p * p;
            if (start < (long) p * p)
                start = (long) p * p;
            for (long x = start; x < highExclusive; x += p) {
                isPrime[(int) (x - low)] = false;
            }
        }

        if (low == 0) {
            isPrime[0] = false;
            if (size > 1)
                isPrime[1] = false;
        } else if (low == 1) {
            isPrime[0] = false;
        }

        List<Long> primes = new ArrayList<>();
        for (int i = 0; i < size; i++) {
            if (isPrime[i])
                primes.add(low + i);
        }
        return primes;
    }

    static long D(long a, long b, long k) {
        long m = k - 1;
        List<Long> primes = segmentedPrimes(a, a + b);
        long total = 0;
        for (long p : primes) {
            total += modInverse(m, p);
        }
        return total;
    }

    public static void main(String[] args) {
        long a = 1000000000L;
        long b = 100000L;
        long k = 100000L;
        System.out.println(D(a, b, k));
    }
}