Problem 779: Prime Factor and Exponent

Project Euler Problem 779 Solution

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

Problem Summary For each integer $n>1$, let $s(n)$ be its smallest prime factor and let $a(n)=v_{s(n)}(n)$ be the exponent of that prime in $n$. The quantity to compute is the limiting average $$A=\lim_{N\to\infty}\frac{1}{N-1}\sum_{n=2}^{N}\frac{a(n)-1}{s(n)-1}.$$ The key observation is that integers can be grouped by their smallest prime factor and by the exact exponent of that prime. That turns the limit into a convergent series over primes, which is what the implementations evaluate. Mathematical Approach We work with natural densities. For a fixed prime $p$, define $$D(p)=\prod_{q<p}\left(1-\frac{1}{q}\right).$$ This is the density of integers not divisible by any smaller prime, so it is the correct prefactor when the smallest prime factor is exactly $p$. Step 1: Density of a Fixed Smallest Prime Factor If $s(n)=p$, then no prime $q<p$ divides $n$. In natural density, each condition $q\nmid n$ contributes a factor $1-1/q$, so the density of integers surviving all smaller primes is $$D(p)=\prod_{q<p}\left(1-\frac{1}{q}\right).$$ For example, $D(2)=1$, $D(3)=1/2$, and $D(5)=(1-1/2)(1-1/3)=1/3$. Step 2: Distribution of the Exponent Now fix $p$ and ask for the density of integers with $s(n)=p$ and $a(n)=k$. Such integers are divisible by $p^k$ but not by $p^{k+1}$....

Detailed mathematical approach

Problem Summary

For each integer $n>1$, let $s(n)$ be its smallest prime factor and let $a(n)=v_{s(n)}(n)$ be the exponent of that prime in $n$. The quantity to compute is the limiting average

$$A=\lim_{N\to\infty}\frac{1}{N-1}\sum_{n=2}^{N}\frac{a(n)-1}{s(n)-1}.$$

The key observation is that integers can be grouped by their smallest prime factor and by the exact exponent of that prime. That turns the limit into a convergent series over primes, which is what the implementations evaluate.

Mathematical Approach

We work with natural densities. For a fixed prime $p$, define

$$D(p)=\prod_{q<p}\left(1-\frac{1}{q}\right).$$

This is the density of integers not divisible by any smaller prime, so it is the correct prefactor when the smallest prime factor is exactly $p$.

Step 1: Density of a Fixed Smallest Prime Factor

If $s(n)=p$, then no prime $q<p$ divides $n$. In natural density, each condition $q\nmid n$ contributes a factor $1-1/q$, so the density of integers surviving all smaller primes is

$$D(p)=\prod_{q<p}\left(1-\frac{1}{q}\right).$$

For example, $D(2)=1$, $D(3)=1/2$, and $D(5)=(1-1/2)(1-1/3)=1/3$.

Step 2: Distribution of the Exponent

Now fix $p$ and ask for the density of integers with $s(n)=p$ and $a(n)=k$. Such integers are divisible by $p^k$ but not by $p^{k+1}$. The $p$-part therefore contributes

$$\frac{1}{p^k}-\frac{1}{p^{k+1}}=\frac{p-1}{p^{k+1}}.$$

Multiplying by the smaller-prime factor gives

$$\delta_{p,k}=D(p)\frac{p-1}{p^{k+1}},\qquad k\ge 1.$$

This is the density of integers whose smallest prime factor is $p$ and whose exponent of $p$ is exactly $k$.

Step 3: Convert the Average into a Prime Series

Each such integer contributes $(k-1)/(p-1)$ to the mean, so

$$A=\sum_{p}\sum_{k\ge 1}\delta_{p,k}\frac{k-1}{p-1}=\sum_{p}D(p)\sum_{k\ge 1}\frac{k-1}{p^{k+1}}.$$

Write $m=k-1$. Then

$$\sum_{k\ge 1}\frac{k-1}{p^{k+1}}=\frac{1}{p^2}\sum_{m\ge 0}\frac{m}{p^m}.$$

Using the standard identity $\sum_{m\ge 0}mx^m=x/(1-x)^2$ with $x=1/p$, we get

$$\sum_{k\ge 1}\frac{k-1}{p^{k+1}}=\frac{1}{p^2}\cdot\frac{1/p}{(1-1/p)^2}=\frac{1}{p(p-1)^2}.$$

Therefore the desired constant is

$$\boxed{A=\sum_{p}\frac{D(p)}{p(p-1)^2}=\sum_{p}\frac{1}{p(p-1)^2}\prod_{q<p}\left(1-\frac{1}{q}\right).}$$

Step 4: Companion Series Used for Verification

The implementations also track the auxiliary mean obtained by replacing the denominator $p-1$ with $p$:

$$A_1=\lim_{N\to\infty}\frac{1}{N-1}\sum_{n=2}^{N}\frac{a(n)-1}{s(n)}.$$

The same density calculation yields

$$A_1=\sum_{p}D(p)\sum_{k\ge 1}\frac{k-1}{p}\cdot\frac{p-1}{p^{k+1}}=\sum_{p}\frac{D(p)}{p^2(p-1)}.$$

This companion series is numerically checked against a direct average over a large initial range of integers, which gives an independent confirmation of the density argument.

Step 5: Worked Example with the First Prime Terms

The first few contributions to $A$ are easy to compute explicitly.

For $p=2$, we have $D(2)=1$, so the contribution is

$$\frac{1}{2(2-1)^2}=\frac{1}{2}.$$

For $p=3$, the prefactor is $D(3)=1/2$, so the contribution is

$$\frac{1/2}{3(3-1)^2}=\frac{1}{24}.$$

For $p=5$, the prefactor is $D(5)=1/3$, so the contribution is

$$\frac{1/3}{5(5-1)^2}=\frac{1}{240}.$$

After only three primes, the partial sum is

$$\frac{1}{2}+\frac{1}{24}+\frac{1}{240}=\frac{131}{240}\approx 0.5458333333,$$

which already shows how quickly the series converges.

Step 6: Why Truncating at a Large Prime Bound Works

Because $0<D(p)\le 1$, each omitted term beyond a cutoff $B$ is at most $1/(p(p-1)^2)$, which behaves like $1/p^3$. So the tail converges absolutely and is extremely small.

The implementation uses the explicit estimate

$$R_B<\frac{1}{2B^2}.$$

With $B=2\cdot 10^6$, this is below $5\times 10^{-13}$, which is comfortably smaller than the displayed $12$-decimal precision.

How the Code Works

The C++, Python, and Java implementations first generate all primes up to $2\cdot 10^6$ with a standard sieve. They then maintain the running product $D(p)$ in increasing prime order, updating it by multiplying with $(p-1)/p$ after each prime has been processed.

At every prime, the implementation adds the target term $D(p)/(p(p-1)^2)$ to the main total. It also accumulates the companion term $D(p)/(p^2(p-1))$, because that auxiliary constant is used as a numerical cross-check of the derivation.

The C++ implementation performs two additional sanity checks. First, it verifies that the explicit tail estimate at the chosen cutoff is below $5\times 10^{-13}$. Second, it checks the companion series both against the known decimal value $0.282419756159$ and against a direct average over the first $200000$ integers obtained by extracting the smallest prime factor and its exponent. The Python and Java implementations keep the same series evaluation but omit the extra assertions and simply print the target value to $12$ decimal places.

Complexity Analysis

Let $B=2\cdot 10^6$ be the sieve bound. Building the prime list costs $O(B\log\log B)$ time and $O(B)$ memory. The subsequent summation visits each prime once, so it costs $O(\pi(B))$ arithmetic steps. High-precision decimal operations increase constant factors but do not change the asymptotic bounds. Overall, the method uses $O(B)$ memory and runs in $O(B\log\log B)$ time.

Footnotes and References

Problem page: Project Euler 779
$p$-adic valuation: Wikipedia — p-adic valuation
Prime factors: Wikipedia — Prime factor
Natural density: Wikipedia — Natural density
Mertens' theorems: Wikipedia — Mertens' theorems
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes

Problem 779 source code

C++

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

#include <boost/multiprecision/cpp_dec_float.hpp>

using boost::multiprecision::cpp_dec_float_100;

static std::vector<int> sieve_primes(int n) {
    std::vector<bool> is_prime(n + 1, true);
    is_prime[0] = false;
    if (n >= 1) {
        is_prime[1] = false;
    }
    for (int i = 2; 1LL * i * i <= n; ++i) {
        if (!is_prime[i]) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            is_prime[j] = false;
        }
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[i]) {
            primes.push_back(i);
        }
    }
    return primes;
}

static long double brute_f1_average(int n) {
    std::vector<int> spf(n + 1, 0);
    for (int i = 2; i <= n; ++i) {
        if (spf[i] != 0) {
            continue;
        }
        spf[i] = i;
        if (1LL * i * i > n) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            if (spf[j] == 0) {
                spf[j] = i;
            }
        }
    }
    long double sum = 0.0L;
    for (int x = 2; x <= n; ++x) {
        int p = spf[x];
        int t = x;
        int a = 0;
        while (t % p == 0) {
            t /= p;
            ++a;
        }
        sum += static_cast<long double>(a - 1) / static_cast<long double>(p);
    }
    return sum / static_cast<long double>(n - 1);
}

int main() {
    using Real = cpp_dec_float_100;

    const int max_prime = 2'000'000;
    const auto primes = sieve_primes(max_prime);

    Real prefix = 1;
    Real f1 = 0;
    Real total = 0;

    for (int p_int : primes) {
        const Real p = Real(p_int);
        const Real pm1 = p - 1;

        f1 += prefix / (p * p * pm1);
        total += prefix / (p * pm1 * pm1);
        prefix *= pm1 / p;
    }

    const Real tail_bound = Real(1) / (Real(2) * Real(max_prime) * Real(max_prime));
    assert(tail_bound < Real("5e-13"));

    const Real stated_f1 = Real("0.282419756159");
    assert(abs(f1 - stated_f1) < Real("1e-9"));

    const long double brute = brute_f1_average(200000);
    const long double theo = f1.convert_to<long double>();
    assert(std::fabsl(brute - theo) < 2e-4L);

    std::cout << std::fixed << std::setprecision(12) << total << '\n';
    return 0;
}

Python

import decimal

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

def solve():
    decimal.getcontext().prec = 50
    primes = sieve_primes(2000000)
    
    prefix = decimal.Decimal(1)
    f1 = decimal.Decimal(0)
    total = decimal.Decimal(0)
    
    for p_int in primes:
        p = decimal.Decimal(p_int)
        pm1 = p - 1
        
        f1 += prefix / (p * p * pm1)
        total += prefix / (p * pm1 * pm1)
        prefix *= pm1 / p
        
    return f"{total:.12f}"

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

Java

import java.math.BigDecimal;
import java.math.MathContext;
import java.util.ArrayList;

public class Euler779 {

    static ArrayList<Integer> sievePrimes(int n) {
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;

        for (int i = 2; i * i <= n; i++) {
            if (isPrime[i]) {
                for (int j = i * i; j <= n; j += i) {
                    isPrime[j] = false;
                }
            }
        }

        ArrayList<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    public static String solve() {
        MathContext mc = new MathContext(50);
        ArrayList<Integer> primes = sievePrimes(2000000);

        BigDecimal prefix = BigDecimal.ONE;
        BigDecimal f1 = BigDecimal.ZERO;
        BigDecimal total = BigDecimal.ZERO;

        for (int pInt : primes) {
            BigDecimal p = new BigDecimal(pInt);
            BigDecimal pm1 = p.subtract(BigDecimal.ONE);

            BigDecimal termF1 = prefix.divide(p.multiply(p).multiply(pm1), mc);
            f1 = f1.add(termF1, mc);

            BigDecimal termTotal = prefix.divide(p.multiply(pm1).multiply(pm1), mc);
            total = total.add(termTotal, mc);

            prefix = prefix.multiply(pm1).divide(p, mc);
        }

        // Output with 12 decimal places
        return String.format(java.util.Locale.US, "%.12f", total);
    }

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