This problem can be solved by a program that finds a prefix sum of a multiplicative function.

What is a multiplicative function?

A function \(f(n)\) that takes a positive integer is said to be a multiplicative function if and only if:

• \(f(ab) = f(a) f(b)\) for all positive integers \(\gcd(a, b) = 1\) with \(\gcd(a, b) = 1\).

In other words, we can interpret that if a positive integer \(n\) can be represented as \(n = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}\) using primes \(p_1, p_2, \dots, p_k\), then \(f(n)\) can be represented as the product of \(f(p_1^{e_1}), f(p_2^{e_2}), \dots, f(p_k^{e_k})\).

A popular problem on a multiplicative function is computation of a prefix sum, specifically:

Prefix sum of multiplicative function

Given a multiplicative function \(f(n)\), evaluate the following expression (assuming \(F(p^e)\) can be evaluated in \(\mathrm{O}(1)\) time for any prime \(p\)):

\[\sum_{1 \leq n \leq N} f(n).\]

The original problem can be actually boiled down to a prefix sum problem of a multiplicative function.

We first assert several facts.

Fact 1

Let \(\sigma(n)\) be the function that returns the sum of the divisors of \(n\). Then \(\sigma(n)\) is multiplicative.

This fact follows from the formula

\[\sigma(n) = \prod_{i=1}^k (1 + p_i + \dots + p_i^{e_i}),\]

where \(n = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}\).

Fact 2

Given a positive integer \(M\), let \(g(n)\) be the function that returns the number of length-\(M\) sequences of positive integers such that the total product of the elements is \(M\). Then \(g(n)\) is multiplicative.

This fact follows from the combinatoric observation showing

\[g(n) = \prod_{i=1}^k \binom{e_i + M - 1}{M - 1},\]

where \(n = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}\). (We leave the proof to the readers; suitable for blue or yellow coders)

Therefore, the original problem can be rephrased as follows:

The original problem reworded

Given a multiplicative function \(g(n)\), evaluate the following expression (assuming \(g(p^e)\) can be evaluated in \(\mathrm{O}(1)\) time for any prime \(p\)):

\[\sum_{1 \leq n \leq N, \sigma(n) \bmod 3 = 0} g(n)\]

If it were not for \(\sigma(n) \bmod 3 = 0\), it is merely a prefix sum problem of a multiplicative function. Can we somehow get rid of it?

Define a multiplicative function \(h(n)\) so that the following is true for any prime \(p\):

\[ h(p^e) = \begin{cases} 0 & \sigma(p^e) \bmod 3 = 0\\ g(p^e) & \text{otherwise}. \\ \end{cases} \]

Then the multiplicativity of \(\sigma(n)\) yields:

\[ h(n) = \begin{cases} 0 & \sigma(n) \bmod 3 = 0 \\ g(n) & \sigma(n) \bmod 3 \neq 0. \\ \end{cases} \]

Therefore, the problem can can be further rephrased as follows:

The original problem reworded (2)

Given a multiplicative function \(g(n)\) and \(h(n)\), evaluate the following expression (assuming \(g(p^e)\) and \(h(p^e)\) can be evaluated in \(\mathrm{O}(1)\) time for any prime \(p\)):

\[\sum_{n=1}^N \left( g(n) - h(n) \right)\]

Of course the expression can be deformed as

\[\sum_{n=1}^N g(n) - \sum_{n=1}^N h(n),\]

which is a prefix sum of \(g(n)\) subtracted by a prefix sum of \(h(n)\). Thus, the original problem has been boiled down to prefix sum problem of multiplicative functions.

There are several ways to compute prefix sum of a multiplicative function. Here we introduce solutions using algorithms called:

• Lucy DP, and
• the black algorithm or Min_25 sieve.

Lucy DP

From now on, we will consider the prefix sum of a multiplicative function \(f(n)\) over \(1 \leq n \leq N\). Define functions \(F(n)\) and \(F_{\mathrm{prime}}(n)\), and integer sets \(q_n\) as follows:

\[F(n) = \sum_{1 \leq m \leq n} f(m),\]

\[F_{\mathrm{prime}}(n) = \sum_{1 \leq p \leq n, p : \mathrm{prime}} f(p),\]

\[Q_n = \left\lbrace m : m = \left \lfloor \frac{n}{k} \right \text{There exists an integer }k\text{ with }\rfloor.\right\rbrace\]

We will try to compute the prefix sum of the multiplicative function by the following two steps:

1. Evaluate \(F_{\mathrm{prime}}(n)\) for all \(n \in Q_N\).
2. Evaluate \(F(N)\) using \(F_{\mathrm{prime}}(n)\).

In the former step, we use a DP (Dynamic Programming) called Lucy DP; in the latter, we use black algorithm or Min_25 sieve.

We first explain Lucy DP. Lucy DP originates from a post by Lucy_Hedgehog in the thread of Project Euler: Problem 10 (accessible after solving).

For simplicity, assume that \(f(p) = 1\) holds for all primes \(p\) as an example. In this case, \(F_{\mathrm{prime}}(n)\) coincides with the prime-counting function \(\pi(n)\), so we will explain how to evaluate \(\pi(n)\).

For \(1 \leq x \leq \lfloor \sqrt{N} \rfloor\) and \(n \in Q_N\), define \(S(x, n)\) as follows:

• \(S(x, n)\): the number of integers between \(2\) and \(n\) that were not sieved by the primes up to \(x\).

Since \(\pi(n) = S(\lfloor \sqrt{N} \rfloor, n)\), it is sufficient to evaluate all \(S(\lfloor \sqrt{N} \rfloor, \ast)\).

Consider representing \(S(x, n)\) using \(S(x-1, \ast)\) to derive the transition of DP.

• If \(x\) is not a prime

No integers will be sieved by \(x\), so \(S(x,n) = S(x-1,n)\).

• If \(x\) is a prime satisfying \(n \lt x^2\):

Multiples of \(x\) up to \(n\) have been already sieved; again, \(S(x, n) = S(x-1, n)\).

• If \(x\) is a prime satisfying \(n \geq x^2\)

The integers \(n\) for which \(nx\) is sieved by \(x\) form the set obtained by excluding the set of primes less than \(x\) from the set of integers not sieved by any primes less than \(x\). Thus,

\[S(x, n) = S(x-1,n) - \left(S(x-1, \lfloor n/x \rfloor) - \pi(x-1) \right);\]

noticing \(\pi(x-1) = S(x-1,x-1)\), we have

\[S(x, n) = S(x-1,n) + S(x-1, \lfloor n/x \rfloor) + S(x-1, x-1).\]

Therefore, the following equation holds:

\[ S(x, n) = \begin{cases} S(x-1, n) + S(x-1, \lfloor n/x \rfloor) + S(x-1, x-1) & (x \text{ is prime}) \wedge x^2 \leq n \\ S(x-1, n) & \mathrm{otherwise} \end{cases} \]

This DP can be updated in-place, with the only indices being \(n\).

• Sample implementation of the function evaluating \(\pi(N)\) in Python

import math
def Pi(N):
    sqrtN = math.isqrt(N)
    Q = [N // i for i in range(1, sqrtN + 1)]
    Q += list(range(Q[-1] - 1, 0, -1))
    S = {i: i - 1 for i in Q}
    for x in range(2, sqrtN + 1):
        if S[x] > S[x - 1]:
            for n in Q:
                if n < x * x: break
                S[n] -= S[n // x] - S[x - 1]
    return S[N]

Now we will estimate the complexity of the DP. We will assume the observations 1, 2, 3 in the editorial of ABC239Ex. For example, \(|Q_N| = \mathrm{O}(\sqrt{N})\).
Now, it is sufficient to count how many times the line S(v) -= S[v // p] - S[p - 1] is executed. We will split cases depending on which of \(x\) and \(N^{\frac{1}{4}}\) is larger.

• If \(2 \leq x \leq N^{\frac{1}{4}}\)

By the Prime Number Theorem, there are \(\mathrm{O}\left(\frac{N^{\frac{1}{4}}}{\log N}\right)\) primes, and the size of the DP table is \(\mathrm{O}(\sqrt{N})\), so it is executed \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) times in total.

• If \(N^{\frac{1}{4}} \leq x \leq \sqrt{x}\)

For each \(x\), we update the parts satisfying \(x^2 \leq n \leq N\). Since \(\sqrt{N} \leq x^2\), the update occurs \(\mathrm{O}\left(\frac{N}{x^2}\right)\) times. Thus, together with the Prime Number Theorem, it is executed \(\mathrm{O}\left(\frac{1}{\log N} \sum_{x=N^{1/4}}^{\sqrt{N}} \frac{N}{x^2} \right) = \mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) times.

Hence, the overall complexity is \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\).

Here, we have explained using \(\pi(n)\), but this algorithm can be used to numerate \(F_{\mathrm{prime}}(n)\) for other multiplicative functions too.
The behavior of the functions \(g\) and \(h\) in our problem against primes \(p\) is as follows:

\[g(p) = M,\]

\[h(p) = \begin{cases} 0 & p \bmod 3 = 2 \\ M & \mathrm{otherwise.} \end{cases} \]

Since \(g(p)\) is constant, the sought value can be found by enumerating \(\pi(n)\) and multiplying them by \(M\). \(h(p)\) is a bit difficult, but they can be basically evaluated by counting primes congruent to \(1\) and \(2\) modulo \(3\) using Lucy DP (details omitted).

The black algorithm

Now we have evaluated \(F_{\mathrm{prime}}(n)\) for \(n \in Q_N\). Now we will explain black algorithm, which evaluates \(F(N)\) using these results.

• As in inventor’s article, it is notable that this algorithm is simple.

First, consider an \(N\)-vertex rooted tree satisfying the following conditions. Here, for \(n \geq 2\), we denote by \(\mathrm{gpf}(n)\) the greatest prime factor of \(n\).

• Its vertices are numbered \(1\) through \(N\).
• Vertex \(1\) is the root, and the parent of vertex \(n\) \((n \geq 2)\) is \(\frac{n}{\mathrm{gpf}(n)}\).

For example, the tree for \(N = 16\) looks like this. (Circles represent leaves, and squares represent non-leaves. The color of the edge connecting vertex \(n\) \((n \geq 2)\) and its parent corresponds to \(\mathrm{gpf}(n)\).)

One property of this tree is that the path from a vertex \(n\) to the vertex \(1\) contains edges corresponding to the prime factors of \(n\). Thus, by performing a DFS (Depth-First Search) appropriately on this tree while maintaining the current \(f(n)\), one can enumerate \(f(1), f(2), \dots, f(N)\) in \(\mathrm{O}(N)\) times, and the prefix sum in \(\mathrm{O}(N)\) too - which is slow.

Let us optimize this DFS algorithm using \(F_{\mathrm{prime}}(n)\). Now, suppose that \(f(n)\) is known for a non-leaf vertex \(n\). If we let

• \(p_k\) be the \(k\)-th smallest prime,
• \(p_i = \mathrm{gpf}(n)\), and
• \(p_j\) be the largest prime not greater than \(N/n\),

then the children of \(n\) are \(n p_i, n p_{i+1}, \dots, n p_j\) in ascending order. Can we find the sum of \(f(n p_i), f(n p_{i+1}), \dots, f(n p_j)\) at once? Putting \(n p_i\) aside as a special case, the primes \(p_{i+1}, \dots, p_j\) are coprime with \(n\), so by the multiplicativity of \(f\) we have

\[\sum_{k=i+1}^j f(n p_k) = f(n) \sum_{k=i+1}^j f(p_k) = f(n) \left( F_{\mathrm{prime}} \left( \left \lfloor \frac{N}{n} \right \rfloor \right) - F_{\mathrm{prime}}(p_i) \right),\]

which enables to instantly find the total sum using \(F_{\mathrm{prime}}(n)\). Using this equation, by finding the sum of \(f\) of children every time entering a vertex \(n\), we do no longer have to visit leaf vertices in the DFS. (For the \(N = 16\) example, we only have to visit vertices \(n=1,2,3,4,8\).)

By appropriately implementing the algorithm above, \(F(N)\) can be evaluated in \(\mathrm{O}(\text{the number of non-leaf vertices in the rooted tree})\). The complexity can be mathematically written as \(\mathrm{O}(N^{1 - \varepsilon})\) (which means: \(\mathrm{O}(N)\), and \(\Omega(N^{1 - \varepsilon})\) for all \(\varepsilon \gt 0\)); but actually, for \(N = 10^{10}\) the number of non-leaves is no more than \(6298637 \simeq 6.3 \times 10^6\), which implies that the constant factor is very lightweight. (The black algorithm works fast enough for the common constraints in competitive programming of \(N \leq 10^{12}\).)

While the black algorithm is a concise algorithm, but meanwhile disadvantageous in that the complexity is \(\mathrm{O}(N^{1 - \varepsilon})\). As an alternative approach, we also introduce the algorithm called Min_25 sieve.

Min_25 sieve

Min_25 sieve is an algorithm proposed by Min_25. Unlike the black algorithm, Min_25 sieve has a complexity evaluation of \(\mathrm{O}(N^{\frac{2}{3}})\), which is \(\mathrm{o}(N)\).
Also, while the black algorithm only evaluates \(F(N)\), Min_25 sieve also yields \(F(n)\) for \(n \in Q_N\) as a byproduct.
Nonetheless, the original Min_25 sieve is a very complex algorithm, so here we introduce a simplified \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) algorithm.

Intuitively, Min_25 sieve does an inverse of Lucy DP. For \(1 \leq x \leq \lfloor \sqrt{N} \rfloor\) and \(n \in Q_N\), define \(S(x, n)\) as follows:

• \(S(x, n)\): he sum of \(f(n)\) over all integers between \(2\) and \(n\), inclusive, that were not sieved by primes not greater than \(x\) in Eratosthenes’ sieve.

Obviously, \(S(\lfloor \sqrt{N} \rfloor, n) = F_{\mathrm{prime}}(n)\), and what we want is \(S(1, N) + 1\). Thus, let us try to compute in the opposite direction of Lucy DP, starting from \(S(\lfloor \sqrt{N} \rfloor, \ast)\) and decreasing the index.

Just as Lucy DP, consider the difference between \(S(x, n)\) and \(S(x-1, n)\).

• If \(x\) is not a prime

Just as Lucy DP, \(S(x, n) = S(x-1, n)\).

• If \(x\) is a prime satisfying \(n \lt x^2\):

Just as Lucy DP, \(S(x, n) = S(x-1, n)\).

• If \(x\) is a prime satisfying \(n \geq x^2\)

Any integers sieved by \(x\) has a form \(x^c m\) (where \(m\) is an integer that can be represented as a product of primes greater than \(x\)). Fixing the value \(c\), the possible \(n\) are

• in the set of the integers not greater than \(n/x^c\) that were remaining when sieved by primes up to \(x\),
• but not in the set of primes up to \(x\),
• or \(1\) if \(c \geq 2\).

Processing the equation based on this fact, we arrive at the following equation:

\[S(x-1, n) = S(x, n) + \sum_{1 \leq c, x^{c+1} \leq n} f(x^c)\left( S(x,\lfloor n/x^c \rfloor) - F_{\mathrm{prime}}(x) \right) + f(x^{c+1}).\]

All that left is to evaluate the recurrence relations based on the equation above. Just as Lucy DP, it can be optimized as an in-place DP. The complexity can also be evaluated in the same manner as Lucy DP, turning out to be \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\).

Both this algorithm and Lucy DP can be optimized using Fenwick tree, reducing the complexity to \(\mathrm{O}\left(N^{\frac{2}{3}}\right)\); the true Min_25 algorithm adopts this optimization. (The author doesn’t know the details though…)

By implementing the algorithms described so far, i.e. Lucy DP and [the black algorithm or Min_25 sieve], the prefix sum of \(g(n)\) and \(h(n)\) can be evaluated, so the problem can be solved. The complexity is \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\), or \(\mathrm{O}(N^{1 - \varepsilon})\) with a small constant factor, which are fast enough.

• [Sample code using the simplified Min_25 sieve (Python)(https://atcoder.jp/contests/abc370/submissions/57434000)]

Advanced topic: even faster algorithm

In this algorithm, we introduced two algorithm that finds a prefix sum of a multiplicative function, the black algorithm and Min_25 sieve; but there are other such algorithms. Investigation is especially active in the Chinese competitive programming community, giving rise to many algorithms.
In terms of complexity, the \(\mathrm{\tilde{O}}(\sqrt{N})\) algorithm developed by zhoukangyang is the fastest. (blog, thesis) The fastest practical algorithm suitable for competitive programming seems to be the \(\Theta \left( \frac{N^{\frac{2}{3}} }{\left(\log N\right)^{\frac{4}{3}}}\right)\) algorithm by Zhenting Zhu.
If the multiplicative function satisfies a certain condition, one can also use Dujiao sieve or PowerfulNumber sieve.
Many people proposes algorithms to find a prefix sum of a multiplicative function, so try searching if you are interested.

References

• Library Checker : Problem proposal (Sum of Multiplicative Function): curated list of past researches on prefix sums of multiplicative functions.
• OI-Wiki: Dujiao sieve and PowerfulNumber sieve are explained.