亚线性筛
时间复杂度可以达到小于线性时间，可以解决 (1e8-1e11) 范围内的问题。比如：杜教筛，(min 25) 筛，州阁筛。
作用
符号规定：(p) 表示素数，(P) 表示素数集合
求大部分积性函数的的前缀和问题，形如：
可以看成一个动态规划的模型，解决很多质因子的贡献问题。
复杂度
(O(frac{n^{frac{3}{4}}}{log n}))
条件
- (f(p)) 可以表示为简单多项式；
- (f(p^k)) 可以快速求；
大致思路
将结果分为两部分：一个是质数部分的和，一个是合数部分的和。
质数部分：
首先，对每一个不同的 (x=lfloor frac{n}{i} floor)，求 (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^2leq n) 时，我们需要考虑哪些数被筛掉了。被筛掉的数的最小质因子一定是 (P_j)，考虑减去 (f(P_j) imes 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=lfloorfrac{n}{i} floor) （此处的 (n) 表示题目中的数据范围）代入到公式中的 (n)，一个 (x) 代表一个块，求多次，目的是为了更快的求出 (g(n,|P|))。
应用：
求 ([1,n]) 素数的个数，(2leq 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
", 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
", g[id]), 0;
}
区间素数和：http://acm.hdu.edu.cn/showproblem.php?pid=6889
求 (sum_{i=3}^{n+1}i+sum_{i=3}^{n+1}i[iin Prime])的值即为结果。
(1leq n leq 10^{10},10^8leq 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
",0);
continue;
}
ans=(n+2)%k*(n+1)%k*inv%k;
ans=(ans+solve()%k-5+k)%k;
printf("%lld
",ans);
}
return 0;
}
合数部分
在素数部分中，我们已经求出了，对于 (x=lfloorfrac{n}{i} floor)，(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)=pigoplus c (p 为质数，igoplus 表示异或))；
- (f(ab)=f(a)f(b) (a,b互质))；
求 (sum_{i=1}^{n}{f(i)} mod (10^9+7))，(1leq n leq 10^{10})。
分析
除 (2) 以外，(f(p)=p-1)。
令 (g(n,j)=sum_{i=1}^{n}{i [iin P|min(i)>P_j]})，那么 (g(n,|P|)=sum_{i=1}^{n}{i[iin P]})。
令 (h(n,j)=sum_{i=1}^{n}{1 [iin P | min(i)>P_j]})，那么 (h(n,|P|)=sum_{i=1}^{n}{1[iin 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
", (S(n, 1) + 1) % mod);
return 0;
}
// 9999998765 986477040
// 9919260817 677875815
定义积性函数 (f(x))，且 (f(p^k)=p^k(p^k-1))，其中 (p) 为素数，求：
对 (1e9+7) 取模。
(1leq 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
",(S(n,1)+1)%mod);
return 0;
}
参考博客：https://www.cnblogs.com/zhoushuyu/p/9187319.html
注意：由于数据比较大，所以很容易爆 ( ext{long long})，有乘法一定要注意取模。