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:
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:
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) \):
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:
\(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];
inline ll read() {
#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() {
n = read(),k = read(),P = read(); Mod.init(P);
init();
for(ll i = 1;i <= sq;i ++) ans = (ans + i * i % Mod * g[i] % Mod) % Mod;
printf("%lld\n",ans);
return 0;
}
Part 2 answer
Or dp. Consider using the first half for reference to define dp array \ (S(n,j) \):
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:
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:
\(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,
Sp1[cnt] = Add(Sp1[cnt - 1],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();
printf("%lld\n",Add(S(n,0),1));
return 0;
}
Make \ (\ sigma_0(n) \) the divisor of \ (n \). Given \ (n,k \), find:
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:
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;
}