min25 筛

亚线性筛

时间复杂度可以达到小于线性时间，可以解决 (1e8-1e11) 范围内的问题。比如：杜教筛，(min 25) 筛，州阁筛。

作用

符号规定：(p) 表示素数，(P) 表示素数集合

求大部分积性函数的的前缀和问题，形如：

[sum_{i=1}^{n}{f(i)} ]

可以看成一个动态规划的模型，解决很多质因子的贡献问题。

复杂度

(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)=sum_{i=1}^{n}{f(i)[iin P or min(p)>P_j,p|i,pin P]} ]

即 (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)=egin{cases} g(n,j-1)&,P_j^2>n\ \ g(n,j-1)-f(P_j) imes[g(frac{n}{P_j},j-1)-sum_{i=1}^{j-1}f(P_i)] &,P_j^2leq n\ end{cases} ]

(g(n,j)) 的初值：(g(n,0)=sum_{i=1}^{n}f(i))，即把所有数都当作是质数带入 (f(p)) 的那个多项式中算出的结果。

因为最后要求得是 (g(n,|P|)) ，因此求得时候数组只需要开第一维，将第二维简化掉。

至此，我们可以利用质数部分来求 ([1,n]) 内的质数的个数，其中 (f(i)=1)。

由此，原公式可以化简为：

[g(n,j)=egin{cases} g(n,j-1)&,P_j^2>n\ \ g(n,j-1)-[g(frac{n}{P_j},j-1)-(j-1)] &,P_j^2leq n\ end{cases} ]

把每个 (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)

将公式变换为：

[g(n,j)=egin{cases} g(n,j-1)&,P_j^2>n\ \ g(n,j-1)-P_j imes[g(frac{n}{P_j},j-1)-sum(j-1)] &,P_j^2leq n\ end{cases} ]

其中，(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)})。

设

[S(n,j)=sum_{i=1}^{n}{[min(p)geq P_j]f(i)} ]

即最小质因子大于等于 (P_j) 的 (f) 值之和。

那么最终的答案为：

[S(n,1)+f(1) ]

因为 (S(n,1)) 中并没有计算 (f(1)) 的值，所以要加上。

因为质数的答案已经算出来了，为 (g(n,j)-sum_{i=1}^{j-1}{f(P_i)})。联系 (g(n,j)) 的含义，要保证最小质因子大于等于 (P_j) ，因此要把小于它的质数的 (f) 值减掉。然后加上合数部分的答案即可。

最终的结果为：

[S(n,j)=g(n,|P|)-sum_{i=1}^{j-1}{f(P_i)+sum_{k=j}^{P_k^2leq n}{sum_{e=1}^{P_k^{e+1}leq n}{S(frac{n}{P^e_{k}},k+1) imes f(P^e_k)+f(P_k^{e+1})}}} ]

其中，(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|]) 的质数进行倍数的枚举，先枚举最小质因子，然后枚举其指数。

模板

loj6053. 简单的函数

定义 (f(x))：

1. (f(1)=1)；
2. (f(p^c)=pigoplus c (p 为质数，igoplus 表示异或))；
3. (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

洛谷PP5325

定义积性函数 (f(x))，且 (f(p^k)=p^k(p^k-1))，其中 (p) 为素数，求：

[sum_{i=1}^{n}f(i) ]

对 (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})，有乘法一定要注意取模。