box's blog

By box, 3 years ago, In English

Consider the following problem (general form of original Min25 sieve):

We are given some multiplicative function $$$f(n)$$$, described at prime powers by $$$f(p^k)=g(p,k)$$$, where $$$g(p,1)$$$ is a polynomial of low degree, and we can calculate $$$g(p,k)$$$ for any singular $$$p$$$ and $$$k$$$ in $$$O(1)$$$.

How do we calculate the prefix sum of $$$f(n)$$$ at the distinct positive integer values of $$$\lfloor\frac Nk\rfloor$$$, where $$$N\le 10^{10}$$$? Clearly, brute force or regular sieves are not going to cut it here.

Prerequisite: For any integer $$$N$$$, for all integers $$$1\le k\le N$$$ there are at most $$$2\sqrt N$$$ distinct values of $$$\lfloor\frac Nk\rfloor$$$. Call these $$$2\sqrt N$$$ such values the blocks of $$$N$$$. (A note on notation: In the following text, $$$\frac Nk$$$ will be used to describe regular division, and $$$N/k$$$ will be used to denote floored division.)


Our workflow is going to be similar to this:

  1. We iterate over each nonzero exponent $$$T$$$ of $$$g(n,1)$$$. We first use the prefix sum of $$$x^T$$$ to estimate the contribution of $$$x^T$$$ to each prefix sum.
  2. Then, we sieve away composite values to ensure that our new prefix sum only counts prime $$$n$$$.
  3. After summing the result of the step 3 sieve over all $$$T$$$, we have obtained the prefix sum of $$$f(n)[n\in\mathbb P]$$$, where $$$\mathbb P$$$ is the set of primes.
  4. Then, we sieve back composite values to restore a correct answer.

The first step is done trivially by an interpolation of your choice.

The second step asks us to evaluate the following function at the blocks of $$$N$$$:

$$$S(n)=\sum\limits_{i=1}^n[i\in\mathbb P]i^T$$$

Let $$$p_1,p_2,\dots$$$ be the sequence of primes, $$$p_0=1$$$, and $$$\mathbb P(i)$$$ be the minimum prime factor of $$$i$$$.

Consider the following, modified, version of $$$S$$$, which "simulates" sieving away primes:

$$$S(n,j)=\sum\limits_{i=1}^n[i\in\mathbb P\lor\mathbb P(i)>p_j]i^T$$$

As $$$j$$$ tends to infinity, each $$$n$$$ satisfies $$$S(n)=S(n,j)$$$. In particular, when $$$p_j^2>n$$$, we have $$$S(n)=S(n,j)$$$, because if $$$x$$$ is composite its minimum prime factor must not be greater than $$$\sqrt x$$$! This means that we don't need to check all primes up to $$$N$$$, but instead only up to $$$\sqrt N$$$. Now, we have the edge case:


The problem becomes: how do we convert $$$S(n,j-1)$$$ evaluated at the blocks of $$$N$$$ to $$$S(n,j)$$$ evaluated at the blocks of $$$N$$$? We effectively need to sieve away the composites with a minimum prime factor of exactly $$$p_j$$$.

In fact, the integers with a minimum prime factor of exactly $$$p_j$$$ are part of $$$p_j^TS(n/p_j,j-1)$$$, because this fixes at least one power of $$$p_j$$$, but unfortunately this includes primes less than $$$p_j$$$. So we simply subtract away them, and get what we need to sieve away: $$$p_j^T(S(n/p_j,j-1)-S(p_{j-1},j-1))$$$.

This gives us an recurrence:


Note that $$$S(p_{j-1},j-1)$$$ is simply the sum of the $$$T$$$-th powers of the first $$$j-1$$$ primes, so this can be precalculated. We replace this with $$$W_{j-1}$$$:


It is now apparent that the values of $$$S(n,j)$$$ at blocks of $$$N$$$ only depend on values of $$$S(n,j-1)$$$ at blocks of $$$N$$$; hence we maintain $$$S(n,j)$$$ at blocks of $$$N$$$ as we slowly increase $$$j$$$ until $$$p_j^2>N$$$.

What is the time complexity of this?

Observe that for some $$$j$$$, we only need to update values of $$$n$$$ where $$$p_j^2<n$$$, because at other positions, $$$S(n,j-1)$$$ is already a correct estimate of $$$S(n)$$$. (Without this optimization, the time complexity becomes $$$O(N/\log N)$$$.)

This means that the time complexity of this part is proportional to the sum of the number of primes less than $$$n$$$ for each $$$n$$$ in the blocks of $$$N$$$. Here it suffices to consider only the $$$n>\sqrt N$$$, as brute force works with the other half.

We have:

$$$\sum\limits_{k=1}^{\sqrt N}\pi(\sqrt{N/k})=\sum\limits_{k=1}^{\sqrt N}O(\frac{\sqrt{\frac Nk}}{\log\frac Nk})$$$
$$$=O(\frac{\int_{1}^{\sqrt N}\sqrt{\frac Nk}dk}{\log N})=O(\frac{N^{3/4}}{\log N})$$$

Now, we move on to the fourth step.

To sieve composite values back, we consider "reversing" the process we did in the second step. This suggests that we still use a function sieving based on the minimum prime factor. Define

$$$R(n,j)=\sum\limits_{i=2}^n[i\in\mathbb P\lor\mathbb P(i)\ge p_j]f(n)$$$

Our base case is now

$$$R(n,\pi(\sqrt n)+1)=S(n)$$$

Our problem becomes to transition from $$$j+1$$$ to $$$j$$$; that is, to sieve back in integers with a minimum prime factor of $$$p_j$$$.

The composite number in question, $$$x$$$, could fall into two categories.

  1. $$$x$$$ may have prime factors that are not $$$p_j$$$; in other words, $$$x=p_j^cy$$$ for some $$$1\le c$$$ and $$$p_j<\mathbb P(y)$$$.
  2. $$$x$$$ is a non-prime power of $$$p_j$$$; in other words, $$$x=p_j^c$$$ for some $$$1<c$$$. Note that, if $$$c=1$$$, then $$$x$$$ is not composite.

We count the contribution of all $$$c$$$ where $$$p_j^c\le n$$$. For the second type, we may count by evaluating $$$g(p_j,c)$$$ with brute force, so we care about deriving the value of the first type.

We must select $$$y$$$ such that $$$1\le y\le n/p_j^c$$$ and $$$\mathbb P(y)>p_j$$$. Like our calculation for the second step, this quantity is $$$R(n/p_j^c,j+1)-\text{included primes}=R(n/p_j^c,j+1)-S(p_j)$$$. If we aren't careful in our enumeration of $$$c$$$, it might be that $$$p_j>n/p_j^c$$$, so we take $$$\max(0,R(n/p_j^c)-S(p_j))$$$.

Remember that $$$f$$$ is multiplicative, so to account for the power $$$p_j^c$$$, we just multiply this by $$$g(p_j,c)$$$.

To summarize what we are sieving in:

$$$\sum\limits_{c=1,p_m^c\le n}g(p_j,c)\max(0,R(n/p_j^c)-S(p_j))+\sum\limits_{c=2,p_m^c\le n}g(p_j,c)$$$

Observe that for the first summation, $$$c$$$ contributes if and only if $$$p_m^{c+1}\le n$$$, so we can rewrite the summation and remove the maximum operand — to get a full recurrence:

$$$R(n,j)=R(n,j+1)+\sum\limits_{c=1,p_m^{c+1}\le n}g(p_j,c)(R(n/p_j^c)-S(p_j))+g(p_j,c+1)$$$

Iterating $$$j$$$ from $$$\pi(\sqrt n)+1$$$ to $$$0$$$, our answers are at $$$R(n,0)$$$.

The analysis of the time complexity of this part is practically the same as that of the first part.


(Yes, in fact this entire blog post was to advertise this problem of mine.)

First let's weaken the problem a bit: How do we calculate

$$$g(n,k,r)=\sum\limits_{i=1}^ni^k[L(i)\equiv r\pmod 4]$$$

where $$$L(i)$$$ is the integer log function, equal to the sum of the prime divisors of $$$i$$$, counting repetition?

At first it is not obvious at all how this is a multiplicative function. However, the Min25 sieve is not bound to functions in the integers. It works perfectly fine for polynomials too!

Consider the multiplicative function $$$\mathfrak L(i)$$$, defined by $$$\mathfrak L(p^c)=p^{ck}x^{(pc)\bmod 4}$$$, operating on $$$\mathbb Z/p\mathbb Z[x]/(x^4-1)$$$. This happens to be the cyclical convolution ring of polynomials of degree $$$4$$$ — in other words, convolution here is defined so that $$$x^{i+j}$$$ overflows to $$$x^{(i+j)\bmod 4}$$$. This is exactly what we need — addition modulo $$$4$$$!

Unfortunately, this function does not seem to satisfy an important requirement of the sieve: $$$p^kx^{p\bmod 4}$$$ is not a polynomial, nor a fully multiplicative function. A first instinct would be to do FFT/NTT — but this does not help either.

The solution would be to completely circumvent the first step. Instead, we take the true meaning of the function $$$p^kx^{p\mod 4}$$$. To find the prefix sum of this at $$$n$$$, what we really want is for each $$$z=0,1,2,3$$$, the sum of the $$$k$$$-th powers of all primes less than $$$n$$$, where said prime modulo $$$4$$$ is equivalent to $$$z$$$.

Surprisingly, this is doable.

To formalize, we need, at the blocks of $$$N$$$:

$$$z=0,1,2,3:\sum\limits_{i=1}^{n}[i\in\mathbb P\land i\equiv z\pmod 4]i^k$$$

Observe that for $$$z=0$$$ the answer is always $$$0$$$, and for $$$z=2$$$ the answer is $$$[2\le n]2^k$$$. Consider how to sieve away the answers for $$$z=1$$$ and $$$z=3$$$ in parallel.

Let us recall the recurrence for the second step of the sieve:


In order to incorporate $$$z$$$, we need to augment this into a polynomial, with nonzero coefficients of $$$x^1$$$ and $$$x^3$$$. Let the current prime in question be $$$p=p_j$$$, and let the previous prime be $$$p'=p_{j-1}$$$. We want

$$$[x^{z|z\in\{1,3\}}]S(n,j)=\sum\limits_{i=1}^n[i\in\mathbb P\lor\mathbb P(i)>p_j][i\equiv z\pmod 4]i^k$$$

Remember that from $$$p'$$$ to $$$p$$$, we sieve away the composites with a minimum prime factor of $$$p$$$, so these composites can all be expressed in the form $$$pd$$$. Now, we can do casework on $$$p$$$ modulo $$$4$$$. In the following, $$$p,pd,d$$$ are all expressed modulo $$$4$$$.

  1. If $$$p=1$$$, then:
    1. For $$$pd=1$$$, we must have $$$d=1$$$. Hence we seive away $$$p^k[x^1](S(n/p_j,j-1)-W_{j-1})$$$ from $$$[x^1]S(n,j)$$$.
    2. For $$$pd=3$$$, we must have $$$d=3$$$. Hence we seive away $$$p^k[x^3](S(n/p_j,j-1)-W_{j-1})$$$ from $$$[x^3]S(n,j)$$$.
  2. If $$$p=3$$$, then:
    1. For $$$pd=1$$$, we must have $$$d=3$$$. Hence we seive away $$$p^k[x^3](S(n/p_j,j-1)-W_{j-1})$$$ from $$$[x^1]S(n,j)$$$.
    2. For $$$pd=3$$$, we must have $$$d=1$$$. Hence we seive away $$$p^k[x^1](S(n/p_j,j-1)-W_{j-1})$$$ from $$$[x^3]S(n,j)$$$.

To calculate $$$S(n,0)$$$, use any polynomial interpolation algorithm.

This recurrence now allows us to successfully complete the second step of sieving. The fourth step is virtually the same, except we operate in $$$\mathbb Z/p\mathbb Z[x]/(x^4-1)$$$.

Notice that the given function at primes has not changed in the original problem. We only need to modify the multiplication operation in the fourth step to complete our solution.

I hope you enjoyed this blog! If you have any interesting problems that can be solved with this sieve, let me know.

Have a nice day :)

  • Vote: I like it
  • +161
  • Vote: I do not like it

| Write comment?
3 years ago, # |
  Vote: I like it 0 Vote: I do not like it

Is this a variation on the derivation of prime-counting (Meissel-Lehmer) algorithm?

  • »
    3 years ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    I guess it counts as an extension. Meissel-Lehmer is a closely related algorithm that doesn't need step 3 and step 4, making it more amenable to optimizations.

3 years ago, # |
  Vote: I like it 0 Vote: I do not like it

It's a standard trick in this type of number theory harmonic-lemma algorithm to use the brute-force technique all the way up to $$$\tilde{O}(n^{2/3})$$$ instead of only up to $$$\sqrt{n}$$$. This results in a slightly better time complexity, $$$\tilde{O}(n^{2/3})$$$ instead of $$$\tilde{O}(n^{3/4})$$$.

2 years ago, # |
  Vote: I like it +46 Vote: I do not like it

Who is this Min25 guy anyway? They're like a modern Prometheus, bringing the most arcane NT-algorithms into our measly CP realm. Yet I've never heard anything about them except their nickname.

2 years ago, # |
  Vote: I like it +13 Vote: I do not like it

The idea of using two different functions that are only equal at prime arguments is brilliant!

Also there is a simpler description of step 2 that shows how to solve this problem for any number instead of just 4.
What we need is to find $$$\sum f(p a_i)$$$ given $$$p$$$ and $$$\sum f(a_i)$$$. In our case, $$$f(n)[x] = n^k x^n$$$ -- a polynomial in $$$x$$$.
But then $$$f(p a)[x] = (p a)^k x^{p a} = p^k (a^k (x^p)^a) = p^k f(a)[x^p]$$$ which is linear at can be applied to the whole sum.

17 months ago, # |
  Vote: I like it 0 Vote: I do not like it

Hi! Not directly related, but I recently noticed that Min25's blog is cleared and empty. The old link is (from their twitter), do you by any chance know what happened?

5 months ago, # |
Rev. 2   Vote: I like it +8 Vote: I do not like it

Sorry for necroposting, but there's an interesting (but very inefficient) way to calculate pi using the Min25 sieve.

Consider the prefix sum of the Euler phi function: $$$\Phi(n) := \sum_{i=1}^n \phi(i)$$$. Since the phi function is multiplicative and is defined by $$$\phi(p^k)=p^{k-1}(p-1)$$$ on prime powers, we can compute this using the Min25 sieve.

On the other hand, it's known that $$$\Phi(n) = \frac{3}{\pi^2} n^2 + O(n \log n)$$$ (see this Wikipedia article). Roughly speaking, this means that $$$\Phi(n)$$$ is approximately $$$\frac{3}{\pi^2} n^2$$$. If we let $$$n$$$ be a power of $$$10$$$, we can use its digits to estimate $$$\frac{3}{\pi^2}$$$.

For example, assuming my implementation is correct, I get $$$\Phi(10^{12})=303963550927059804025910$$$. So $$$\frac{3}{\pi^2} \approx 0.303963550927059804025910$$$, from which we get $$$\pi \approx \sqrt{\frac{3}{0.303963550927059804025910}} \approx 3.1415926535895533$$$.