# [learning notes] Min_25 sieve

Immortal min_ The magic sieve invented by the invention is used to sieve the prefix and of the product function.

For the integrable function \ (f(x) \) to sieve the prefix sum, the specific requirements are that \ (f(p) \) is a simple polynomial, and \ (f(p^e) \) can be calculated quickly.

Complexity is sub linear, but I won't prove it.

## Description

For the prefix sum of an integral function \ (f(x) \), we can make the following transformation:

\begin{aligned} \sum_{i=1}^nf(x)=&\sum_{p\le n}f(p)+\sum_{\substack{i\le n\\i\not\in\mathbb P}}f(i)\\ =&\sum_{p\le n}f(p)+\sum_{\substack{p\le \sqrt n\\ p^e\le n}}\sum_{\substack{i\le\lfloor\frac n{p^e}\rfloor\\ \operatorname{lpf}(i)>p}}f(i) \end{aligned}

The first line is to divide the contribution of the number of \ (1\sim n \) into prime numbers and composite numbers for calculation respectively.

The second line is to enumerate the minimum prime factor of each composite number, and then calculate the rest, where \ (\ operatorname{lpf}(i) \) represents the minimum prime factor of \ (i \).

Then it can be divided into two parts to calculate!

### Part 1 first half

Given that \ (f(p) \) is a simple polynomial about \ (P \), we can divide this polynomial into several monomials for calculation. In other words, we only need to find \ (\ sum {P \ Le n} p^k \), where \ (p^k \) is a part of the original polynomial.

Consider a magical dp. Define the dp array \ (g(n,j) \) as:

$g(n,j)=\sum_{i\le n}[i\in \mathbb P\operatorname{or}\operatorname{lpf}(i)>p_j]i^k$

If the algorithm process is imagined as an Ehrlich sieve, \ (g(n,j) \) actually calculates the contribution of all the "primes" remaining after screening \ (p_, J \) and multiples to the answer.

Consider moving from \ (g(n,j-1) \), that is, we need to calculate the contribution of screening out \ (p_, J \) to the answer. Divide \ (p_j \) into two categories:

• $$p_j > \ sqrt n, p_j^2 > n$$: it can be seen that the minimum composite number with the minimum prime factor of \ (p_j \) is \ (p_j^2 \), so screening out \ (p_j \) does not contribute to the answer, and we can get \ (g(n,j)=g(n,j-1) \).
• $$p_j\le\sqrt n,p_j^2\le n$$: enumerate the values of the composite numbers screened out after removing the minimum prime factor \ (p_j \). It can be seen that the contribution of screening out \ (p_j \) to the answer is \ (p_j^kg(\lfloor\frac n{p_j}\rfloor,j-1) \). But in this way, we will screen \ (p_1\sim p_{j-1} \) again, so the answer should add \ (p_j^kg(p_j-1,j-1) \).

So we get the recurrence formula of \ (g(n,j) \):

$g(n,j)=g(n,j-1)-p_j^k \begin{cases} 0&p_j>\sqrt n,p_j^2>n\\ g(\lfloor\frac n{p_j}\rfloor,j-1)-g(p_j-1,j-1)&p_j\le\sqrt n,p_j^2\le n \end{cases}$

Where \ (g(p_j-1,j-1) \) is actually \ (\ sum {I = 1} ^ {J-1} p_j ^ k \), which can be linearly screened out and recorded as \ (\ textit{Sp}_j \).

Note that when \ (p_j > \ sqrt n \), the values of \ (g(n,j) \) are the same, so the value of such \ (g(n,j) \) can be omitted. At the same time, for each \ (x \), we only need the values of all \ (g(\lfloor\frac nd\rfloor,x) \) and only \ (O(\sqrt n) \) such values. We can map the subscript during storage.

Example: P5493 prime prefix statistics

Let \ (S(n) \) be the sum of \ (k \) powers of prime numbers within \ (n \).
Given \ (N \), find:

$\sum_{i=1}^{\lfloor\sqrt N\rfloor}i^2S(\lfloor\frac Ni\rfloor)\bmod p$

$$0 \le k \le 10$$，$$1 \le N \le 4\times {10}^{10}$$，$${10}^9 < p < 1.01 \times {10}^9$$ .

The board problem can be realized according to the above process. Calculate the power of natural number \ (k \) and refer to CF622F.

code (using FastMod to optimize modulus):

#include <cstdio>
#include <cmath>

const int maxn = 1e6 + 5;
typedef long long ll;

namespace FastModnmsp {
typedef unsigned int uint;
typedef __uint128_t uL;
typedef unsigned long long ull;

struct FastMod {
ull b,m;
FastMod() {}
FastMod(ull _b): b(_b),m(ull((uL(1) << 64) / _b)) {}
inline void init(ull _b) { b = _b,m = ull((uL(1) << 64) / _b); }
friend inline ull operator % (const ull& a,const FastMod& mod) {
ull r = a - (uL(mod.m) * a >> 64) * mod.b;
return r >= mod.b ? r - mod.b : r;
}
} Mod;
} using FastModnmsp::Mod;

ll n,sq,k,ans,P,finv[15],prsk[15];
ll cntP,pri[maxn],prik[maxn],sumk[maxn];
ll cntS,w[maxn],g[maxn];
int id1[maxn],id2[maxn];
bool isp[maxn];

#define gc c = getchar()
ll d = 0; int f = 0,gc;
for(;c < 48 || c > 57;gc) f |= (c == '-');
for(;c > 47 && c < 58;gc) d = (d << 1) + (d << 3) + (c ^ 48);
#undef gc
return f ? -d : d;
}

inline ll fpow(ll a,ll b) {
ll res = 1;
for(;b;a = a * a % Mod,b >>= 1) if(b & 1) res = res * a % Mod;
return res;
}

inline ll getPowSum(ll n) {
ll ans = 0,p = 1,pre = 0;
for(int i = 1;i <= k + 2;i ++) {
p = p * (n - i) % Mod,pre = (pre + prsk[i]) % Mod;
if(n == i) return pre;
ll inv = finv[i - 1] * finv[k + 2 - i] % Mod * fpow(n - i,P - 2) % Mod;
ll res = pre * inv % Mod;
if((k - i) & 1) res = P - res;
ans = (ans + res) % Mod;
}
return p * ans % Mod;
}

inline int getId(ll k) { return k <= sq ? id1[k] : id2[n / k]; }

inline void init() {
finv[0] = isp[1] = 1; sq = sqrt(n);
for(ll i = 1;i <= k + 2;i ++) finv[i] = finv[i - 1] * fpow(i,P - 2) % Mod,prsk[i] = fpow(i,k);
for(ll i = 2;i <= sq;i ++) {
if(!isp[i]) pri[++ cntP] = i,
prik[cntP] = fpow(i,k),
sumk[cntP] = (sumk[cntP - 1] + prik[cntP]) % Mod;
for(int j = 1;j <= cntP && i * pri[j] <= sq;j ++) {
isp[i * pri[j]] = true;
if(!(i % pri[j])) break;
}
}
for(ll r,l = 1;l <= n;l = r + 1) {
r = n / (n / l); w[++ cntS] = n / l;
g[cntS] = getPowSum(w[cntS] % Mod) - 1;
if(n / l <= sq) id1[n / l] = cntS;
else id2[r] = cntS;
}
for(int i = 1;i <= cntP;i ++)
for(int j = 1;j <= cntS && pri[i] * pri[i] <= w[j];j ++) {
ll h = getId(w[j] / pri[i]);
g[j] -= prik[i] * (g[h] - sumk[i - 1] + P) % Mod;
if(g[j] < 0) g[j] += P;
}
}

int main() {
init();
for(ll i = 1;i <= sq;i ++) ans = (ans + i * i % Mod * g[i] % Mod) % Mod;
printf("%lld\n",ans);
return 0;
}

Or dp. Consider using the first half for reference to define dp array \ (S(n,j) \):

$S(n,j)=\sum_{i\le n}[\operatorname{lpf}(i)>p_j]f(i)$

The answer is obviously \ (S(n,0)+1 \). Here \ (1 \) is the contribution of \ (f(1) \), which is proposed for calculation.

Next, the contribution is discussed in two parts:

• Prime: the contribution of all prime numbers is \ (g(n,j) \), but the contribution of \ (p_1\sim p_{j-1} \) should be subtracted, so the answer is \ (g(n,j)-\textit{Sp}_j\).
• Composite number: enumerate the minimum prime factor and its number of times to obtain the contribution of \ (\ sum {\ substack {k > J \ \ p_k^e \ Le n}} f (p_k^e) (s (\ lfloor \ frac n{p_k^e}, K) + [e\not=1]) \), where \ ([e\not=1] \) calculates the contribution of \ (f(1) \) after removing \ (p_k^e \).

So the final formula is:

$S(n,j)=g(n,x)-\textit{Sp}_j+\sum_{\substack{k>j\\p_k^e\le n}}f(p_k^e)(S(\lfloor\frac n{p_k^e},k)+[e\not=1])$

Where \ (x \) is the largest \ (x \) within \ (n \), that is \ (x=\pi(\lfloor\sqrt n\rfloor) \).

It's amazing that this thing can be searched directly. It doesn't need to be memorized at all, and the complexity is excellent \ (O (\ frac {n ^ {0.75} {\ log n}) \).

Example:

Let the integrable function \ (f(p^k)=p^k(p^k-1) \). Given \ (n \), find:

$\sum_{i=1}^nf(i)\bmod (10^9+7)$

$$1\le n\le 10^{10}$$.

At the prime number, there is obviously \ (f(p)=p(p-1)=p^2-p \), which can be Min_25 sieve.

code:

#include <cstdio>
#include <cmath>

typedef long long ll;
const int maxn = 1e6 + 5;
const ll P = 1e9 + 7,inv3 = (P + 1) / 3;

ll n,sq,tot,cnt;
ll pri[maxn]; bool isp[maxn];
ll Sp1[maxn],Sp2[maxn];
ll g1[maxn],g2[maxn],w[maxn];
int id1[maxn],id2[maxn];

inline ll getLiSum(ll n) { return n * (n + 1) / 2 % P; }
inline ll getSqSum(ll n) { return n * (n + 1) / 2 % P * (2 * n + 1) % P * inv3 % P; }

inline ll Add(ll a,ll b) { a += b; return a >= P ? a - P : a; }

inline int getId(ll k) { return k <= sq ? id1[k] : id2[n / k]; }

inline void Init() {
isp[1] = true,sq = sqrt(n);
for(ll i = 2;i <= sq;i ++) {
if(!isp[i]) pri[++ cnt] = i,
Sp2[cnt] = Add(Sp2[cnt - 1],i * i % P);
for(int j = 1;j <= cnt && i * pri[j] <= sq;j ++) {
isp[i * pri[j]] = true;
if(!(i % pri[j])) break;
}
}
for(ll r,l = 1;l <= n;l = r + 1) {
r = n / (n / l);
w[++ tot] = n / l;
g1[tot] = getLiSum(w[tot] % P) - 1;
g2[tot] = getSqSum(w[tot] % P) - 1;
if(n / l <= sq) id1[n / l] = tot;
else id2[r] = tot;
}
for(int i = 1;i <= cnt;i ++)
for(int j = 1;j <= tot && pri[i] <= w[j] / pri[i];j ++) {
int k = getId(w[j] / pri[i]);
g1[j] -= pri[i] * (g1[k] - Sp1[i - 1] + P) % P;
g2[j] -= pri[i] * pri[i] % P * (g2[k] - Sp2[i - 1] + P) % P;
if(g1[j] <= 0) g1[j] += P; if(g2[j] <= 0) g2[j] += P;
}
}

inline ll S(ll n,int j) {
if(pri[j] >= n) return 0;
ll x = getId(n),ans = Add(((g2[x] - Sp2[j]) - (g1[x] - Sp1[j])) % P,P);
for(int k = j + 1;k <= cnt && pri[k] <= n / pri[k];k ++) {
ll p = pri[k];
for(int e = 1;p <= n;e ++,p *= pri[k]) {
ll tmp = p % P;
ans = Add(ans,tmp * (tmp - 1) % P * (S(n / p,k) + (e != 1)) % P);
}
}
return ans;
}

int main() {
scanf("%lld",&n); Init();
return 0;
}

Make \ (\ sigma_0(n) \) the divisor of \ (n \). Given \ (n,k \), find:

$\sum_{i=1}^n\sigma_0(i^k)\bmod 2^{64}$

Multi test, \ (1\le n,k\le10^{10} \).

Three times the experience problem. After writing this problem, you can cut DIVCNT2 and DIVCNT3.

Make \ (f (n) = \ Sigma_ 0 (n ^ k \), then obviously:

$f(n)=\begin{cases} 1&n=1\\ ke+1&n=p^k\\ f(a)f(b)&n=ab,a\perp b \end{cases}$

Yes, Min_25 sift it out.

code:

#include <cstring>
#include <cstdio>
#include <cmath>

typedef unsigned long long ull;
const int maxn = 1e6 + 5;

ull n,K,sq;
ull pri[maxn]; bool isp[maxn];
ull Sp[maxn],g[maxn],w[maxn];
int t,tot,cnt,id1[maxn],id2[maxn];

inline int getId(ull k) { return k <= sq ? id1[k] : id2[n / k]; }

inline void Init() {
isp[1] = true,sq = sqrt(n);
for(ull i = 2;i <= sq;i ++) {
if(!isp[i]) pri[++ cnt] = i,
Sp[cnt] = Sp[cnt - 1] + K + 1;
for(int j = 1;j <= cnt && i * pri[j] <= sq;j ++) {
isp[i * pri[j]] = true;
if(!(i % pri[j])) break;
}
}
for(ull r,l = 1;l <= n;l = r + 1) {
r = n / (n / l);
w[++ tot] = n / l;
g[tot] = w[tot] - 1;
if(n / l <= sq) id1[n / l] = tot;
else id2[r] = tot;
}
for(int i = 1;i <= cnt;i ++)
for(int j = 1;j <= tot && pri[i] <= w[j] / pri[i];j ++) {
int k = getId(w[j] / pri[i]);
g[j] -= g[k] - i + 1;
}
for(int i = 1;i <= tot;i ++) g[i] *= (K + 1);
}

inline ull S(ull n,int j) {
if(pri[j] >= n) return 0;
ull x = getId(n),ans = g[x] - Sp[j];
for(int k = j + 1;k <= cnt && pri[k] <= n / pri[k];k ++) {
ull p = pri[k];
for(ull e = 1;p <= n;e ++,p *= pri[k])
ans += (K * e + 1) * (S(n / p,k) + (e != 1));
}
return ans;
}

int main() {
scanf("%d",&t);
for(;t;t --) {
memset(isp,false,sizeof(isp)); tot = cnt = 0;
scanf("%llu%llu",&n,&K); Init();
printf("%llu\n",S(n,0) + 1);
}
return 0;
}

Keywords: Math OI

Added by sheckel on Fri, 21 Jan 2022 00:43:04 +0200