min25 筛
亚线性筛
时间复杂度可以达到小于线性时间，可以解决 \(1e8-1e11\) 范围内的问题。比如：杜教筛，\(\min 25\) 筛，州阁筛。
作用
符号规定：\(p\) 表示素数，\(P\) 表示素数集合
求大部分积性函数的的前缀和问题，形如：
可以看成一个动态规划的模型，解决很多质因子的贡献问题。
复杂度
\(O(\frac{n^{\frac{3}{4}}}{log\ n})\)
条件
- \(f(p)\) 可以表示为简单多项式；
- \(f(p^k)\) 可以快速求；
大致思路
将结果分为两部分：一个是质数部分的和，一个是合数部分的和。
质数部分：
首先，对每一个不同的 \(x=\lfloor \frac{n}{i} \rfloor\)，求 \(\sum_{i=1}^{x}{f(i)}\) [\(i\) 是质数]。先筛出 \(\sqrt{n}\) 以内的质数集合 \(P\)，\(P_j\) 表示从小到大第 \(j\) 个质数。
设：
即 \(g(n,j)\) 表示从\(1\) 到 \(n\) 的 \(f(i)\) 的累加和，其中 \(i\) 为素数或者 \(i\) 的最小质因子大于第 \(j\) 个素数。
考虑该式子的实际意义，帮助理解：
假设 \(f(i)=id(i)\) 即 \(f(i)=i\)，根据埃式筛的原理：每次筛到一个素数，都把它的倍数筛掉。那么 ，\(g(n,j)\) 就表示用埃式筛将第 \(j\) 个素数及其倍数筛出后剩余的所有数之和加上前 \(j\) 个素数的和。
显然 \(ans=g(n,|P|)\)，\(|P|\) 为素数集合的大小，下面考虑如何去求 \(g(n,j)\)，分两种情况：
- 当 \(P_j^2>n\) 时，满足最小质因子为 \(P_j\) 的最小合数为 \(P_j^2 >n\) ，此时第 \(j\) 次筛不会再筛掉任何数，所以：\(g(n,j)=g(n,j-1)\)。
- 当 \(P_j^2\leq n\) 时，我们需要考虑哪些数被筛掉了。被筛掉的数的最小质因子一定是 \(P_j\)，考虑减去 \(f(P_j)\times g(\frac{n}{P_j},j-1)\) （积性函数），但是可以发现 \(g(\frac{n}{P_j},j-1)\) 多减去了那些最小质因子小于 \(P_j\) 的数的 \(f\) 值，因此要加上 \(\sum_{i=1}^{j-1}f(P_i)\)。
综上所述，可得转移方程为：
\(g(n,j)\) 的初值：\(g(n,0)=\sum_{i=1}^{n}f(i)\)，即把所有数都当作是质数带入 \(f(p)\) 的那个多项式中算出的结果。
因为最后要求得是 \(g(n,|P|)\) ，因此求得时候数组只需要开第一维，将第二维简化掉。
至此，我们可以利用质数部分来求 \([1,n]\) 内的质数的个数，其中 \(f(i)=1\)。
由此，原公式可以化简为：
把每个 \(x=\lfloor\frac{n}{i} \rfloor\) （此处的 \(n\) 表示题目中的数据范围）代入到公式中的 \(n\)，一个 \(x\) 代表一个块，求多次，目的是为了更快的求出 \(g(n,|P|)\)。
应用：
求 \([1,n]\) 素数的个数，\(2\leq n \leq 10^{11}\)。
\(\min25\) 筛求区间素数个数模板：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9 + 7;
const int N = 1e6 + 5;
int prime[N], tol;
bool vis[N];
ll w[N], g[N], n;
int sqr, id1[N], id2[N]; // id1[],id2[]存x对应的索引
int id(ll x) { return x <= sqr ? id1[x] : id2[n / x]; } //求出是第几个x
void init() { //欧拉筛
sqr = (int)(sqrt(1.0 * n) + 0.5);
tol = 0;
for (int i = 2; i <= sqr; i++) {
if (!vis[i])
prime[++tol] = i;
for (int j = 1; j <= tol && prime[j] * i <= sqr; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0)
break;
}
}
}
ll solve() {
int cnt = 0; //记录块的个数
for (ll i = 1, j = 0; i <= n; i = j + 1) { //分块
j = n / (n / i); //块的右端点
w[++cnt] = n / i; //有cnt个不同的 x=n/i
if (w[cnt] <= sqr) //因为x的值可能很大，所以分开来保存序号
id1[w[cnt]] = cnt;
else
id2[n / w[cnt]] = cnt;
g[cnt] = (w[cnt] - 1); //每一个x对应一个g[]，对于[1,x]内质数的个数，赋初值为 x-1
}
//按照公式进行转移
for (int i = 1; i <= tol; i++) //枚举素数集合
{
for (int j = 1; j <= cnt; j++) //遍历所有的 x
{
if (1LL * prime[i] * prime[i] > w[j])
break;
int k = id(w[j] / prime[i]);
g[j] -= g[k] - (i - 1);
}
}
return g[id(n)];
}
int main() {
scanf("%lld", &n);
init(); //先筛出sqrt(n)内的素数
printf("%lld\n", solve());
return 0;
}
还有另一种写法（借用他人代码）：
#include <cstdio>
#include <math.h>
#define ll long long
const int N = 316300;
ll n, g[N << 1], a[N << 1];
int id, cnt, sn, prime[N];
inline int Id(ll x) { return x <= sn ? x : id - n / x + 1; }
int main() {
scanf("%lld", &n), sn = sqrt(n);
for (ll i = 1; i <= n; i = a[id] + 1) a[++id] = n / (n / i), g[id] = a[id] - 1;
for (int i = 2; i <= sn; ++i)
if (g[i] != g[i - 1]) {
// 这里 i 必然是质数，因为 g[] 是前缀质数个数
// 当 <i 的质数的倍数都被筛去，让 g[] 发生改变的位置只能是下一个质数
// 别忘了 i<=sn 时，ID(i) 就是 i。
prime[++cnt] = i;
ll sq = (ll)i * i;
for (int j = id; a[j] >= sq; --j) g[j] -= g[Id(a[j] / i)] - (cnt - 1);
}
return printf("%lld\n", g[id]), 0;
}
区间素数和：http://acm.hdu.edu.cn/showproblem.php?pid=6889
求 \(\sum_{i=3}^{n+1}i+\sum_{i=3}^{n+1}i[i\in Prime]\)的值即为结果。
\(1\leq n \leq 10^{10},10^8\leq k \leq 10^9\)
将公式变换为：
其中，\(sum(j-1)\) 表示前 \(j-1\) 个质数的前缀和。
注意取模，取模次数太多会被超时。
区间素数和模板：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2e5+5;
ll n,k,inv,w[N],g[N],sum[N];
int sqr,tol,prime[N],id1[N],id2[N];
bool vis[N];
int getid(ll x)
{
return x<=sqr?id1[x]:id2[n/x];
}
ll power(ll a,ll b)
{
ll res=1;
a%=k;
while(b)
{
if(b&1) res=res*a%k;
a=a*a%k;
b>>=1;
}
return res;
}
void init()
{
tol=0;
sum[0]=0;
for(int i=2;i<=sqr;i++)
{
if(!vis[i])
{
prime[++tol]=i;
sum[tol]=(sum[tol-1]+i);
}
for(int j=1;j<=tol&&prime[j]*i<=sqr;j++)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0) break;
}
}
}
ll solve()
{
n++;
sqr=(int)(sqrt(1.0*n)+0.5);
init();
int cnt=0;
for(ll l=1,r=0;l<=n;l=r+1)
{
r=n/(n/l);
w[++cnt]=n/l;
if(w[cnt]<=sqr) id1[w[cnt]]=cnt;
else id2[n/w[cnt]]=cnt;
g[cnt]=(w[cnt]*(w[cnt]+1)/2-1);
}
for(int i=1;i<=tol;i++)
{
for(int j=1;j<=cnt;j++)
{
if(1LL*prime[i]*prime[i]>w[j]) break;
int tk=getid(w[j]/prime[i]);
g[j]=(g[j]-1LL*prime[i]*(g[tk]-sum[i-1]));
}
}
return g[getid(n)];
}
int main()
{
int t;
scanf("%d",&t);
while(t--)
{
scanf("%lld%lld",&n,&k);
inv=power(2LL,k-2);
ll ans=0;
if(n==1)
{
printf("%d\n",0);
continue;
}
ans=(n+2)%k*(n+1)%k*inv%k;
ans=(ans+solve()%k-5+k)%k;
printf("%lld\n",ans);
}
return 0;
}
合数部分
在素数部分中，我们已经求出了，对于 \(x=\lfloor\frac{n}{i} \rfloor\)，\(\sum_{i=1}^{x}{[i 是质数] f(i)}\)。
设
即最小质因子大于等于 \(P_j\) 的 \(f\) 值之和。
那么最终的答案为：
因为 \(S(n,1)\) 中并没有计算 \(f(1)\) 的值，所以要加上。
因为质数的答案已经算出来了，为 \(g(n,j)-\sum_{i=1}^{j-1}{f(P_i)}\)。联系 \(g(n,j)\) 的含义，要保证最小质因子大于等于 \(P_j\) ，因此要把小于它的质数的 \(f\) 值减掉。然后加上合数部分的答案即可。
最终的结果为：
其中，\(g(n,|P|)\) 表示 \(\sum_{i=1}^{|P|}{f(P_i)}\) ,减去 \(\sum_{i=1}^{j-1}{f(P_i)}\) 后表示 \(\sum_{i=j}^{|P|}{f(P_i)}\) ，根据 \(S(n,j)\) 的定义，我们还需要求出 \([1,n]\) 内，最小质因子大于 \(P_j\) 的合数的 \(f\) 函数值的和。可以通过对 \(P\) 内的第 \([j,|P|]\) 的质数进行倍数的枚举，先枚举最小质因子，然后枚举其指数。
模板
定义 \(f(x)\)：
- \(f(1)=1\)；
- \(f(p^c)=p\bigoplus c\ (p 为质数，\bigoplus 表示异或)\)；
- \(f(ab)=f(a)f(b)\ (a,b互质)\)；
求 \(\sum_{i=1}^{n}{f(i)} \bmod (10^9+7)\)，\(1\leq n \leq 10^{10}\)。
分析
除 \(2\) 以外，\(f(p)=p-1\)。
令 \(g(n,j)=\sum_{i=1}^{n}{i [i\in P|\min(i)>P_j]}\)，那么 \(g(n,|P|)=\sum_{i=1}^{n}{i[i\in P]}\)。
令 \(h(n,j)=\sum_{i=1}^{n}{1 [i\in P | \min(i)>P_j]}\)，那么 \(h(n,|P|)=\sum_{i=1}^{n}{1[i\in P]}\)。
先将 \(2\) 和其它数一样处理，当 \(j=1\) 时，加上 \(2\)。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9 + 7;
const int N = 1e6 + 5;
const ll inv2 = 500000004;
int prime[N], tol, id1[N], id2[N];
bool vis[N];
int sqr;
ll g[N], w[N], sum[N], h[N], n; // h[]计质数的个数
int getid(ll x) { return x <= sqr ? id1[x] : id2[n / x]; }
void init() {
sqr = (int)(sqrt(n) + 0.5);
tol = 0;
sum[0] = 0;
for (int i = 2; i <= sqr; i++) {
if (!vis[i]) {
prime[++tol] = i;
sum[tol] = (sum[tol - 1] + i) % mod;
}
for (int j = 1; j <= tol && prime[j] * i <= sqr; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0)
break;
}
}
}
void solve() {
init();
int cnt = 0;
for (ll l = 1, r = 0; l <= n; l = r + 1) {
r = n / (n / l);
w[++cnt] = n / l;
if (w[cnt] <= sqr)
id1[w[cnt]] = cnt;
else
id2[n / w[cnt]] = cnt;
g[cnt] = ((w[cnt] + 2) % mod) * ((w[cnt] - 1) % mod) % mod * inv2 % mod;
// w[cnt]要先取模在再乘，否则会爆long long
h[cnt] = (w[cnt] - 1) % mod;
}
for (int i = 1; i <= tol; i++) {
for (int j = 1; j <= cnt && 1LL * prime[i] * prime[i] <= w[j]; j++) {
int k = getid(w[j] / prime[i]);
g[j] = (g[j] - 1LL * prime[i] * (g[k] - sum[i - 1]) % mod + mod) % mod;
h[j] = (h[j] - (h[k] - (i - 1)) + mod) % mod;
}
}
}
ll S(ll x, int y) {
if (x <= 1 || prime[y] > x)
return 0;
int k = getid(x);
ll res = (g[k] - h[k] - sum[y - 1] + y - 1 + mod) % mod;
if (y == 1)
res += 2;
for (int i = y; i <= tol && 1LL * prime[i] * prime[i] <= x; i++) {
ll p1 = prime[i], p2 = 1LL * prime[i] * prime[i];
for (int j = 1; p2 <= x; j++, p1 = p2, p2 *= prime[i])
res = (res + S(x / p1, i + 1) * (prime[i] ^ j) % mod + (prime[i] ^ (j + 1)) % mod) % mod;
}
return res;
}
int main() {
scanf("%lld", &n);
solve();
printf("%lld\n", (S(n, 1) + 1) % mod);
return 0;
}
// 9999998765 986477040
// 9919260817 677875815
定义积性函数 \(f(x)\)，且 \(f(p^k)=p^k(p^k-1)\)，其中 \(p\) 为素数，求：
对 \(1e9+7\) 取模。
\(1\leq n \leq 10^{10}\)
分析
\(f(p)=p^2-p\)，将 \(p^2\) 和 \(p\) 分开来维护，其它的套板子即可。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
const int N=1e6+5;
const int maxn=1e4+100;
const int inv6=166666668;
const int inv2=500000004;
int prime[N],tol,sqr,id1[N],id2[N];
bool vis[N];
ll n,sum1[N],sum2[N],w[N],g1[N],g2[N];//g[N];
int getid(ll x)
{
return x<=sqr?id1[x]:id2[n/x];
}
void init()
{
sqr=(int)(sqrt(n)+0.5);
tol=0;
sum1[0]=sum2[0]=0;
for(int i=2;i<=sqr;i++)
{
if(!vis[i])
{
prime[++tol]=i;
sum1[tol]=(sum1[tol-1]+1LL*i*i%mod)%mod;
sum2[tol]=(sum2[tol-1]+i)%mod;
}
for(int j=1;j<=tol&&prime[j]*i<=sqr;j++)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0) break;
}
}
}
void solve()
{
init();
int cnt=0;
for(ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
w[++cnt]=n/l;
if(w[cnt]<=sqr) id1[w[cnt]]=cnt;
else id2[n/w[cnt]]=cnt;
//g[cnt]=(w[cnt]+1)%mod*(w[cnt]%mod)%mod*((2*w[cnt]+1)%mod)%mod;
//g[cnt]=(g[cnt]*inv6)%mod;
//g[cnt]=(g[cnt]-(w[cnt]+1)%mod*(w[cnt]%mod)%mod*inv2%mod+mod)%mod;
g1[cnt]=((w[cnt]+1)%mod*(w[cnt]%mod)%mod*((2*w[cnt]+1)%mod)%mod*inv6%mod-1+mod)%mod;
//特别注意取模
g2[cnt]=((w[cnt]+1)%mod*(w[cnt]%mod)%mod*inv2%mod-1+mod)%mod;
}
for(int i=1;i<=tol;i++)
{
for(int j=1;j<=cnt&&1LL*prime[i]*prime[i]<=w[j];j++)
{
int k=getid(w[j]/prime[i]);
//g[j]=(g[j]-1LL*prime[i]*(prime[i]-1)%mod*(g[k]-sum[i-1])%mod+mod)%mod;
g1[j]=(g1[j]-1LL*prime[i]*prime[i]%mod*(g1[k]-sum1[i-1])%mod+mod)%mod;
g2[j]=(g2[j]-1LL*prime[i]*(g2[k]-sum2[i-1])%mod+mod)%mod;
}
}
}
ll S(ll x,int y)
{
if(x<=1||prime[y]>x) return 0;
int k=getid(x);
ll res=(g1[k]-g2[k]-sum1[y-1]+sum2[y-1]+mod)%mod;
for(int i=y;i<=tol&&1LL*prime[i]*prime[i]<=x;i++)
{
ll p1=prime[i],p2=1LL*prime[i]*prime[i];
for(int j=1;p2<=x;j++,p1=p2,p2*=prime[i])
res=(res+S(x/p1,i+1)*p1%mod*(p1-1)+p2%mod*((p2-1)%mod)%mod)%mod;
}
return res;
}
int main()
{
scanf("%lld",&n);
solve();
printf("%lld\n",(S(n,1)+1)%mod);
return 0;
}
参考博客：https://www.cnblogs.com/zhoushuyu/p/9187319.html
注意：由于数据比较大，所以很容易爆 \(\text{long long}\)，有乘法一定要注意取模。
浙公网安备 33010602011771号