[study notes] min_25 sieve

Posted by stenk on Wed, 26 Jan 2022 10:46:28 +0100

Some agreement

  1. \(p_k \) represents the prime number of \ (K \), especially \ (p_0=1 \);
  2. Set \ (n/k \) to \ (\ left \lfloor \frac{n}{k}\right \rfloor \);
  3. Let \ (\ text{lpf}(x) \) be the minimum prime factor of the number \ (x \).

What are you doing?

Prefix and sum of integrability function \ (F(i) \).

Push persimmon

Make \ (H (I, J) = \ sum_ {x = 2} ^ I [\ text {LPF} (x) > p_j] f (x) \), then the answer is \ (h(n,0)+F(1)=h(n,0)+1 \) (the first term of the integral function must be \ (1 \), \ (1 \) and any number coprime). The general idea of solution seems to be recursion. You might as well expand it according to the definition first (note that this expansion only includes the composite number that can be decomposed into multiple prime numbers):

\[h(i,j)=\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k) \]

In fact, it is to enumerate the minimum prime factor and its power, and then recurse (it is easy to get that they are coprime).

In addition, there are prime numbers and the power of a single prime number, which is also easier to write:

\[h(i,j)=\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k)+F(p_k^{e+1})\ \ \ \ (p_{k}^{e+1}\le i) \]

Note the following restrictions. It can be found that when \ (p_k ^ {e + 1} > I \), the \ (F(p_k^e)\cdot h(i/p_{k}^e,k) \) must also be zero, so this restriction is compact and reasonable.

Let \ (P(n) \) be the prefix and sum of \ (F \) values of prime numbers less than or equal to \ (n \), then there will be (note that prime numbers only need to enumerate to the \ (\ sqrt n \) level):

\[h(i,j)=P(i)-P(p_j)+\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k)+F(p_k^{e+1})\ \ \ \ (p_{k}^{e+1}\le i) \]

So we just calculate \ (P(n) \)

First of all, we need that \ (F(n) \) is a polynomial of small degree, so that each term disassembled can be regarded as a completely integrated function - more generally, let \ (f_k(n)=n^k \), just calculate the \ (P_k(n) \) of this function, multiply it by the coefficients and add it to get \ (P(n) \)

Let \ (g(i,j) \) be the sum of \ (f \) values of numbers between \ ([2,i] \) and \ (\ text {LPF} > p_j \). In particular, this lump always contains \ (f \) values of prime numbers less than or equal to \ (i \).

Consider recursively finding \ (g(i,j) \) - we can understand \ (g(i,j) \) as the value of \ (f \) excluding the composite number with the minimum prime factor of \ (p_, J \) on the basis of \ (g(i,j-1) \):

\[g(i,j)=g(i,j-1)-f_k(p_j)\cdot\left (g(i/p_j,j-1)-\sum_{q=1}^{j-1}f_k(p_q)\right)\ \ \ \ \ (i\ge p_j^2) \]

Since this is a completely integrable function, there is \ (f(r\cdot p_j)=f(r)\cdot f(p_j) \), and the latter one is equivalent to calculating the legal \ (r \) It can be found that the property of complete integrable function eliminates the operation of enumerating the power of prime numbers.

Note the following conditions \ (i\ge p_j^2 \). In fact, if \ (I < p_j ^ 2 \), it is impossible to find a legal \ (r \), so you can directly \ (g(i,j)\leftarrow g(i,j-1) \).

It is easy to find that \ (P_k(n)=g(n,cnt) \) (\ (cnt \) is the number of prime numbers screened out).

About complexity: \ (\ text{It is said that} \) the complexity of calculating \ (h(i,j) \) is \ (\ mathcal O(n^{1-\epsilon}) \), and the complexity of calculating \ (g(i,j) \) is \ (\ mathcal o \ left (\ frac {n ^ {\ frac {3} {4}} {\ log n} \ right) \) Prove no, it may keep cooing.

code implementation

Refer to a valley Template question.

In the implementation process, the number of \ (i,j \) parameters of the function is less than \ (\ sqrt, n \)\ The reason for (j \) will not be explained. For \ (i \), it can be observed that it is a downward nested structure similar to \ (\ left \lfloor \frac{\left \lfloor \frac{n}{a}\right \rfloor}{b}\right \rfloor \), which can prove that its value belongs to the value set of \ (\ left \lfloor \frac{n}{a}\right \rfloor \), which can not only prove that its number is within \ (\ sqrt n \), It can also be divided into blocks for preprocessing. The following is a perceptual proof: it is advisable to replace the denominator in \ (\ left \lfloor \frac{n}{a}\right \rfloor \) with \ (a '\) so that \ (a'\mid n \) and \ (\ frac{n}{a'}=\left \lfloor \frac{n}{a}\right \rfloor \), so that nesting can be replaced with \ (\ left \lfloor \frac{n}{a'b}\right \rfloor \), and the conclusion is proved.

Therefore, we save the reasonable \ (i \). It should be noted that the value of \ (i \) may be very large, so use \ (i \) to find the part of \ (i\le \sqrt n \) and \ (n/i \) to find the other part (this thing actually has \ (i=n/l \) when dividing blocks, and the corresponding \ (r \) is \ (n/i \)).

Returning to this question, we can see that \ (F(n)=n^2-n \), which can be divided into \ (f_1,f_2 \) and then combined.

#include <cstdio>
#define print(x,y) write(x),putchar(y)

template <class T>
inline T read(const T sample) {
	T x=0; char s; bool f=0;
	while((s=getchar())>'9' || s<'0')
		f |= (s=='-');
	while(s>='0' && s<='9')
		x = (x<<1)+(x<<3)+(s^48),
		s = getchar();
	return f?-x:x;
}

template <class T>
inline void write(T x) {
	static int writ[50],tp=0;
	if(x<0) putchar('-'),x=-x;
	do writ[++tp] = x-x/10*10, x/=10; while(x);
	while(tp) putchar(writ[tp--]^48);
}

#include <cmath>
#include <iostream>
using namespace std;
typedef long long ll;

const int maxn = 1e6+5;
const int mod = 1e9+7;
const int inv = 166666668;

bool is[maxn];
ll n,id[maxn];
int S,f[2][maxn],pc,p[maxn];
int g[2][maxn],tot,ind[2][maxn];

inline int inc(int x,int y) {
	return x+y>=mod?x+y-mod:x+y;
}

inline int dec(int x,int y) {
	return x-y<0?x-y+mod:x-y;
}

void sieve(int n) {
	for(int i=2;i<=n;++i) {
		if(!is[i]) {
			p[++pc] = i;
			f[0][pc] = inc(f[0][pc-1],i);
			f[1][pc] = inc(f[1][pc-1],1ll*i*i%mod);
		}
		for(int j=1; j<=pc && i*p[j]<=n; ++j) {
			is[i*p[j]] = 1;
			if(i%p[j]==0) break;
		}
	}
}

inline int getSum(ll n,bool d) {
	int ret; n %= mod; // qwq
	if(!d) ret = (n*(n+1)>>1)%mod;
	else ret = n*(n+1)%mod*(n*2+1)%mod*inv%mod;
	return (ret-1<0?ret-1+mod:ret-1);
}

void getf() {
	ll l,r,x;
	for(l=1;l<=n;l=r+1) {
		r = min(n,n/(n/l)); 
		id[++tot] = (x=n/l);
		if(x<=S) ind[0][x]=tot;
		else ind[1][n/x]=tot; // for better restoration
		g[0][tot] = getSum(x,0);
		g[1][tot] = getSum(x,1);
	}
    // Recursive g(i,j) with rolling array
	for(int i=1;i<=pc;++i) {
		for(int j=1; j<=tot && 1ll*p[i]*p[i]<=id[j]; ++j) {
			int pre = (id[j]/p[i]<=S)?
					ind[0][id[j]/p[i]]:ind[1][n/(id[j]/p[i])];
			g[0][j] = dec(g[0][j],1ll*p[i]*dec(g[0][pre],f[0][i-1])%mod);
			g[1][j] = dec(g[1][j],1ll*p[i]*p[i]%mod*dec(g[1][pre],f[1][i-1])%mod);
		}
	}
}

int h(ll x,int y) {
	if(x<=p[y]) return 0;
	int ID = (x<=S)?ind[0][x]:ind[1][n/x];
	int ans = dec(dec(g[1][ID],g[0][ID]),dec(f[1][y],f[0][y]));
	for(int i=y+1; i<=pc && 1ll*p[i]*p[i]<=x; ++i) {
		for(ll tmp=p[i];;tmp=tmp*p[i]) {
			if(tmp*p[i]>x) break;
			ans = inc(ans,inc(
				tmp*(tmp-1)%mod*h(x/tmp,i)%mod,
				tmp*p[i]%mod*((tmp*p[i]-1)%mod)%mod
			));
		}
	}
	return ans;
}

int main() {
	n=read(9ll);
	sieve(S=sqrt(n)); getf(); 
	print(inc(h(n,0),1),'\n');
	return 0;
}

Topic summary

However, the goose has only done one now