**hidden**

## TEES - Tutorial: Extended Eratosthenes Sieve

### Extended Eratosthenes Sieve:

### A General Approach to Find the Prefix Sum of Multiplicative Functions

This topic has been discussed about 2-3 years ago in Chinese NOI Winter Camp, and finally they ended in an $O(n^{3/4} / \log n)$ general approach. In practise it's slow due to large constant factor. In the following I will give a general approach with about $O(n^{0.7})$ time and $O(\sqrt n)$ memory, which is far faster than the previous method when $N \le 10^{12}$.

Most of the idea is original authored by Hiroaki Yamanouchi. I learned this from his code.

In the following discussion we use $\sqrt N$ to denote the largest integer $W$ where $W^{2} \le N$; $p, p_{1}, p_{2}, ...$ are all primes.

The parameter $magic \ge \sqrt N$ can be chosen freely. But in (my) practise I always choose $magic = \sqrt N$ (to decrease the memory usage and make the code cache-friendly).

Consider the general question: A multiplicative function $f(n)$, defined as:

$f(n) = \begin{cases} 1 &\text{if } n=1 \\ g(p,e) &\text{if } n=p^{e}, e > 0 \\ f(x)*f(y) &\text{if } n = x * y, gcd(x, y)=1 \end{cases}$

Where $f(p) = g(p, 1)$ is any polynomial of $p$. (Other $g(p, i)$ doesn't affect the problem much.)

Given a positive integer $N$. Calculate the prefix sum of this function: $S = \sum_{1\le i\le N} f(i)$ .

To solve this problem, consider the prime factorization of any integer $M \leq N$:

$M = \prod_{1\le i\le k} p_{i}^{e_{i}}$

$p_{1} \lt p_{2} \lt ... \lt p_{k}, e_{i} \gt 0$

The main idea is based on the following observation:

$M' = M / p_{k}$ will not have any prime factor $\gt magic$.

So, we can use following pseudo-code to calculate $S$:

$S = f(1)$ for each $M' \le N$ without any prime factor $\gt magic$: $S += f(M') * \sum_{F \lt p \le N / M'} f(p)$, where $F$ is the largest prime factor of $M'$ and $p$ is a prime if (the exponent of $F$ in $M'$) $\gt 1$: $S += f(M')$ end if end for return S |

Enumerating each $M'$ and calculating $f(M')$ can be easily solved by a simple recursion function. $\sum_{F \lt p \le N / M'} f(p)$ can be viewed as $S'(N/M') - S'(F)$ where $S'(X) = \sum_{1 \lt p \le X} f(X)$ where $p$ is a prime. Now the rest of the problem is: how to calculate $S'(X)$ quickly?

Let's take a look at all the possible values of querying $X$. Since $M$ has no prime factor $\gt magic$, $F$ can only take values of primes $\le magic$. Also, if $M' \gt N / magic$, $N / M' \le magic$. So $N / M'$ can only take values $\le magic$ or $N / i$ where $i \le N / magic$. So we only have to calculate at most $magic + N / magic$ values: $1, 2, 3, ..., ,magic, N / 1, N / 2, ... N / (N / magic)$. We name this set $NS$.

The following pseudo code solves the problem when $f(p)=p^{d}$ where d is any non-negative integer. i.e. Given $N$ and $d$, for each $X$ in set $NS$, calculate $S'(X) = \sum_{1 \lt p \le X} p^{d}$, where $p$ is a prime.

For each $i$ in set $NS$ Map[$i$] = $\sum_{2 \le j \le i} j^{d}$ // This can be calculated by a formula in O(1) End for For $p$ = 2,3,5,7... (any prime not exceed $\sqrt N$) in increasing order: For each $i$ in set $NS$ in decreasing order: If $i \ge p*p$: Map[$i$] -= (Map[$i/p$] - Map[$p-1$]) $*p^{d}$ else break End if End for End for |

This code is indeed simulation of the process of original Eratosthenes Sieve: in each step, we kick out the multiples of prime $p$ which have not been kicked out by any previous primes yet.

When implementing Extended Eratosthenes Sieve we need to store Map[1], Map[2], ..., Map[$magic$] and Map[$N/1$], Map[$N/2$], ... Map[$N / (N/magic)$] into two separated arrays and maintain the array subscripts carefully. When implementing the recursion function to enumerate $M'$, we won't take this kind of $M'$ into consideration:

Suppose the largest prime factor of $M'$ is $F$: $M' * F \ge N$ and $M' \text{ mod } F^{2} \ne 0$.

Note: These two optimizations may affect the analysis of time complexity of the whole algorithm.

Since any polynomial can be written as the weighted sum of several $p^{d}$, to calculate any any polynomial f function with degree d, we can pre-calculate $O(d * magic)$ values. Then when asking for any $S'(X)$, we simply take out the corresponding values and calculate their weighted summation. This completes the whole algorithm.

The memory complexity of this algorithm is obviously $O(magic)$. The time complexity of this algorithm has not been proved by me yet. (According to the performance it approximates function $N^{0.7}$).

The following problem on SPOJ can be used as practise: APS2,DIVCNT3,DIVCNTK.

**Further discussions**

Note 1: According to discussion in Codeforces, for the latter part of the algorithm, we can obtain time complexity about $O(n^{2/3})$. But in practise ($ N \le 10^{12} $) the algorithm above (use $magic = \sqrt N $, with complexity at least $O(n^{3/4} / \log n)$) is faster because of low memory usage (fitting into cache).

Note 2: To obtain time complexity about $O(n^{2/3})$, it seems that the latter part of the algorithm can't be directly calculated. For more information refer to this.

Note 3: The number of $M' \le N$ without any prime factor $\gt \sqrt N$ is about $O(N^{\log 2})$. Please refer to this.

Added by: | Fudan University Problem Setters |

Date: | 2018-01-19 |

Time limit: | 1s |

Source limit: | 50000B |

Memory limit: | 1536MB |

Cluster: | Cube (Intel G860) |

Languages: |