Multiplicative Function and Smallest Prime Factor

When solving Project Euler problems, we often need to compute values such as the number or sum of divisors for every positive integer up to N. Since these functions are easy to evaluate once we know the prime factorization, a natural starting point is to extend the Sieve of Eratosthenes to store the factorization of every number in the range.

let n: usize = 1_000_000;
let mut is_composite: Vec<bool> = vec![false; n];
let mut factors: Vec<Vec<(usize, usize)>> = vec![vec![]; n];

for p in 2..n {
    if is_composite[p]{
       continue;
    }
    // p is prime
    factors[p].push((p, 1));
    for i in (2 * p..=n).step_by(p){
        is_composite[i] = true;

        // save factors
        factors[i].push((p, factor(i, p)));
    }
}

let mut result: Vec<usize> = vec![0; n];
for i in 2..n {
    result[i] = compute(&factors[i]);
}

This approach can produce the answers, but it has two major sources of overhead.

1. Storage: keeping a separate Vec for every number takes substantial space. Around a million heap allocations and scattered memory accesses also add considerable overhead.
2. Repeated computation: we build each factorization and evaluate compute separately. For example, when computing f(999983 × 2), we do not reuse the value of f(999983) that we have already calculated.

We can address the first issue by storing only the Smallest Prime Factor (SPF) and reconstructing the remaining factors when needed. For the second issue, combining SPF with the multiplicative property of the target function lets us reuse previous results through dynamic programming. The Eratosthenes-based SPF sieve still takes O(N log log N) time, but one more refinement gives us an algorithm with linear time complexity. This is called the Linear Sieve. In this article, we implement these approaches and compare their measured performance.

Smallest Prime Factor

The main source of storage overhead in the naive implementation is keeping every prime factor explicitly. This requires an array for each number and, in turn, heap allocations. In fact, knowing just one prime factor is enough to recover the full factorization efficiently. Every integer $n > 1$ has a prime factor $p$, so we can write $n = i p$. If $i != 1$, we can apply the same reasoning to a prime factor $q$ of $i$ and write $n = j p q$, repeating the process as needed. Since the remaining values keep decreasing in the order $n > i > j$, storing one factor per number lets us reuse this information to recover factorizations.

Instead of storing every prime factor during the sieve, we therefore store only the first one found. This is the Smallest Prime Factor approach, implemented below.

let n: usize = 1_000_000;
let mut spf: Vec<usize> = vec![0; n];

for p in 2..n {
    // p is composite
    if spf[p] != 0 {
       continue;
    }
    // p is prime
    for i in (p..n).step_by(p){
        if spf[i] == 0 {
            spf[i] = p
        }
    }
}

fn factor(m: usize, spf: &Vec<usize>) -> Vec<usize> {
    let mut x = m;
    let mut factor_list = Vec::new();
    while x > 1 {
        factor_list.push(spf[x]);
        x /= spf[x];
    }
    factor_list
}

Multiplicative function

An arithmetic function $f: \mathbb{Z} \rightarrow \mathbb{C}$ is multiplicative if it satisfies the following property.

\begin{equation} {}^{\forall}m, n \in \mathbb{Z} \text{ are coprime, then } f(mn) = f(m)f(n) \end{equation}

This property lets us compute the function for every integer using only its values at prime powers.

\begin{equation} \text{if } n = \prod_{i=1}^m p_i^{a_i}, \quad \text{then } f(n) = \prod_{i=1}^m f(p_i^{a_i}) \end{equation}

The divisor count, divisor sum, and Möbius function $\mu$ are examples of multiplicative functions. We can compute their values with the following algorithm.

fn multiplicative_values(
    n: usize,
    value_for_prime_power: impl Fn(usize, usize) -> i64,
) -> Vec<i64> {
    let mut spf = vec![0; n];

    for p in 2..n {
        if spf[p] != 0 {
            continue;
        }
        for m in (p..n).step_by(p) {
            if spf[m] == 0 {
                spf[m] = p;
            }
        }
    }

    let mut f = vec![0; n];
    if n > 1 {
        f[1] = 1;
    }

    for m in 2..n {
        let p = spf[m];
        let (mut rest, mut exponent) = (m, 0);

        while rest % p == 0 {
            rest /= p;
            exponent += 1;
        }

        // m = p^exponent * rest, The two factors are coprime.
        // Since rest < m, f[rest] has already been computed.
        f[m] = value_for_prime_power(p, exponent) * f[rest];
    }

    f
}

By using the smallest prime factor to write $m = p^a r$, we obtain $r < m$. We can therefore compute the value at $m$ using the previously stored function value at $r$.

Linear Sieve

The Sieve of Eratosthenes is already effective, but there is still room for improvement. The loop for m in (p..n).step_by(p) visits every multiple of each prime, so some numbers are visited repeatedly. More precisely, each number is visited once for each of its distinct prime factors. The average number of distinct prime factors grows on the scale of $\log\log N$, giving this sieve a time complexity of $O(N\log\log N)$.

If we could visit each number just once, we could achieve $O(N)$ time. This is the idea behind the linear sieve. Its key observation comes from the definition of SPF. For a number $n$ whose SPF is $p$, write $n = p i$. For $p$ to be the smallest prime factor of $n$, it must be less than or equal to the smallest prime factor of $i$. In other words, $i$ determines the range of primes $p$ that can be used in this decomposition.

So far, we have discovered a prime $p$ and visited its multiples to fill in the sieve. In the linear sieve, we instead visit products of $i$ with only the primes allowed by this condition. The resulting implementation is shown below.

fn linear_sieve(n: usize) -> Vec<usize> {
    let mut spf = vec![0; n];
    let mut primes = Vec::new();

    for i in 2..n {
        if spf[i] == 0 {
            spf[i] = i;
            primes.push(i);
        }

        for &p in &primes {
            // Check the product bound without risking multiplication overflow.
            if p > (n - 1) / i {
                break;
            }

            spf[i * p] = p;

            if p == spf[i] {
                break;
            }
        }
    }

    spf
}

The primes list is sorted in increasing order, so we process only primes satisfying p <= spf[i]. Since the smallest prime factor is unique, the corresponding pair $(i, p)$ for $n$ is uniquely determined as $(n / \operatorname{spf}(n), \operatorname{spf}(n))$. Each composite number is therefore generated exactly once, eliminating repeated visits and giving an $O(N)$ time complexity.

SPF also tells us whether $i$ already contains the factor $p$. If $p < \operatorname{spf}(i)$, then $i$ does not contain $p$, so the two are coprime and we can write $f(n) = f(p)f(i)$. If $p = \operatorname{spf}(i)$, then $i$ is also divisible by $p$. In this case, we need to remove the remaining factors of $p$ from $i$ to obtain coprime factors. By maintaining an array of the exponent of $p$ in each $i$, we can determine how many factors remain and the exponent of $p$ in $n$ without another loop.

The implementation is as follows.

fn multiplicative_values(
    n: usize,
    value_for_prime_power: impl Fn(usize, usize) -> i64,
) -> Vec<i64> {
    let mut spf = vec![0; n];
    let mut primes = Vec::new();

    let mut exponent = vec![0; n];
    let mut rest = vec![0; n];
    let mut memory = vec![0; n];

    memory[1] = 1;
    for i in 2..n {
        if spf[i] == 0 {
            spf[i] = i;
            primes.push(i);

            exponent[i] = 1;
            rest[i] = 1;
            memory[i] = value_for_prime_power(i, 1);
        }

        for &p in &primes {
            // Check the product bound without risking multiplication overflow.
            if p > (n - 1) / i {
                break;
            }

            let x = i * p;
            spf[x] = p;

            if p < spf[i] {
                exponent[x] = 1;
                rest[x] = i;
                memory[x] = memory[p] * memory[i];
            } else if p == spf[i] {
                exponent[x] = exponent[i] + 1;
                rest[x] = rest[i];
                memory[x] = value_for_prime_power(p, exponent[x]) * memory[rest[x]];
                break;
            }
        }
    }
    memory
}

Comparison

To compare the approaches, we compute the divisor count $\tau(m)$ for every integer $1 \le m < N$. Its prime-power rule is $\tau(p^a) = a + 1$, which takes constant time.

The four implementations perform the same task:

• Stored factorizations: build a vector of prime factors and exponents for every number, then multiply the exponent-plus-one terms.
• SPF + factorization: build the Eratosthenes-based SPF table, reconstruct each number’s full factor list as in the SPF example, and count its divisors.
• SPF + DP: use the multiplicative-function implementation to reuse previously computed values.
• Linear sieve + DP: compute the SPF, exponent, remainder, and function-value arrays together.

The measurements below were taken on September 8, 2026, on an AMD Ryzen 7 PRO 8840U, using Rust 1.94.0 on x86-64 Linux and compiling with rustc -O. Each implementation uses one thread. Times are the median of seven runs, with execution order rotated between runs and a correctness/warm-up pass before measurement. The entire size sweep was measured again independently; the table reports the second set of seven-run medians, rather than choosing the faster result from either set.

N (exclusive upper bound)	Stored factorizations	SPF + factorization	SPF + DP	Linear sieve + DP
1,000	0.026 ms	0.026 ms	0.007 ms	0.003 ms
10,000	0.355 ms	0.342 ms	0.094 ms	0.034 ms
100,000	4.196 ms	3.665 ms	0.945 ms	0.353 ms
1,000,000	117.936 ms	40.544 ms	9.963 ms	12.103 ms
10,000,000	1,471.795 ms	540.470 ms	204.843 ms	147.744 ms
100,000,000	19,736.816 ms	7,955.241 ms	2,499.455 ms	1,492.147 ms

Timing includes array allocation, sieve construction, function evaluation, and destruction of temporary storage. It excludes compilation, validation, printing, and destruction of the returned result vector. The benchmark uses black_box to keep the computed results observable. All four implementations were checked against direct divisor counting through 10,000, including small boundary cases, and their complete output arrays were compared at every measured N.

The executable benchmark normalizes the range to 1..N, supplies the divisor-count helpers omitted from the introductory example, changes that example’s inner range from 2 * p..=n to 2 * p..n, and guards the linear implementation’s initialization of memory[1] for small inputs. The snippets above retain their existing executable code; these timings are for the corrected executable benchmark.

At N = 100,000,000, SPF + DP is about 7.9 times faster than storing all factorizations, while linear sieve + DP is about 13.2 times faster. The linear implementation is about 1.68 times faster than SPF

• DP at this size.

The two independent seven-run medians at N = 100,000,000 differed by less than 1% for every method. Within the second set, the minimum and maximum times were:

Method	Minimum	Median	Maximum
Stored factorizations	19.679 s	19.737 s	20.585 s
SPF + factorization	7.719 s	7.955 s	8.062 s
SPF + DP	2.475 s	2.499 s	2.520 s
Linear sieve + DP	1.480 s	1.492 s	1.506 s

Using medians prevents a single unusually fast run from determining the reported result. Repetition does not eliminate every source of measurement bias, but the independent runs support the same performance ordering at this size.

Interestingly, the linear sieve does not win at every size. Both full sweeps showed SPF + DP ahead at N = 1,000,000, despite linear sieve + DP winning at smaller and larger sizes. A finer sweep, also using seven-run medians, shows that this is a region rather than an isolated point:

N	SPF + DP	Linear sieve + DP
200,000	1.926 ms	0.812 ms
400,000	4.103 ms	2.683 ms
600,000	6.018 ms	4.683 ms
800,000	7.825 ms	8.565 ms
1,500,000	16.616 ms	19.517 ms
2,000,000	28.877 ms	25.508 ms
4,000,000	69.579 ms	54.512 ms

In this sweep, linear sieve + DP loses its lead between 600,000 and 800,000, and regains it between 1,500,000 and 2,000,000. The sampled points locate the changes in ordering, not the exact crossover thresholds.

Cache capacity is a plausible explanation. This CPU has a 16 MiB L3 cache. At N = 1,000,000, the two main SPF + DP arrays occupy about 15.3 MiB, while the four main linear-sieve arrays occupy about 30.5 MiB, plus the prime list. The linear implementation therefore reaches this memory scale earlier. At smaller sizes, its lower operation count can dominate; in the middle region, its larger memory footprint may offset that advantage; at larger sizes, avoiding repeated sieve visits can become more valuable again.

This is a hypothesis, not a measured attribution: total array size is not the same as the actively accessed working set, the cache is shared, and the timings also include allocation and temporary-storage destruction. Cache-miss profiling or a controlled change in array representation would be needed to establish the cause. These results describe this implementation and machine and do not imply a universal crossover point.