【数论】Min_25筛

Min_25筛

学习资料：OI Wiki

​ P5325题解 by wucstdio

用法：

求解 \(\sum_{i=1}^n f(i)\) ，满足 \(f(x)\) 是一个积性函数，且 \(f(p^e)\) 是关于 \(p\) 的低阶多项式。

板子&&例题

例题1： P5325 【模板】Min_25筛

题意：

定义积性函数 \(f(x)\) ，且 \(f(p^k)=p^k(p^k-1)\) （ \(p\) 是一个质数)。对给定整数 \(n\) ， \(1\le n\le10^{10}\) ，求 \(\sum_{i=1}^nf(i)\) 。对 \(10^9+7\) 取模。

解：

首先，对 \(i\) 按照质数和合数分类：

\[\sum^n_{i=1}f(i)=\sum_{1\le p\le n}f(p)+\sum^n_{i=1}f(i)[i\notin prime] \]

然后，枚举后边合数的最小质因子以及最小质因子次数。所有合数的最小质因子一定小于等于 \(\sqrt{n}\) ：

\[\sum_{1\le p\le n}f(p)+\sum_{1\le p^e\le n,1\le p,\le{\sqrt{n}}}f(p^e)(\sum_{i=1}^{\frac n{p^e}}f(i)*[minp>p]) \]

其中， \(minp\) 表示 \(i\) 的最小质因子。

将式子分为前后两部分。

\(\clubsuit\) 第一部分：

设完全积性函数 \(g(x)\) ，其在质数处取值 \(g(x)=f(x\))

设数组 \(G(n,j)\) 满足：

\[G(n,j)=\sum^n_{i=1} g(i)[i\in prime\vee minp>p_j] \]

因为 \(G(n,j)\) 可以由 \(G(n,j-1)\) 减去 最小质因子为 \(p_j\) 的合数函数和 得到，而最小质因子为 \(p_j\) 的合数函数和，提取出质因子 \(p_j\) 后可表示为 ：

\[g(j)[G(\frac n{p_j},j)-G(p_{j-1},j-1)] \]

所以有：

\[G(n,j)=G(n,j-1)-g(j)[(G(\frac n{p_j},j-1)-G(p_{j-1},j-1))] \]

其中：

\[G(p_{j-1},j-1)=\sum_{i=1}^{j-1}g(p_j)=\sum_{i=1}^{j-1}f(p_j) \]

所以：

\[G(n,j)=G(n,j-1)-g(p_j)[(G(\frac n{p_j},j-1)-\sum_{i=1}^{j-1}g(p_i)]\\ \]

对于 \(G(n,j)\) 函数，\(1\) 到 \(n\) 中所有质数的函数值和 等值于 \(G(n,x)\) ，其中 \(p_x\) 是最后一个小于等于 \(\sqrt{n}\) 的质数。为了方便，计其为 \(G(n)\) 。可用整除分块求。

\(\clubsuit\) 第二部分：

设 \(S(n,x)\) 表示求 \(1\) 到 \(n\) 中所有最小质因子大于 \(p_x\) 的函数值之和。则题目的答案就是 \(S(n,0)\) （\(p_0=1\) ) 。

将 \(S(n,x)\) 分成两部分，第一部分是大于 \(p_x\) 的质数，即 $G(n)-\sum_{i=1}^x g(p_i) $，第二部分是最小质因子大于 \(p_x\) 的合数，枚举最小质因子：

\[S(n,x)=G(n)-\sum_{i=1}^x g(p_i)+\sum_{p_k^e\le n\wedge k>x}f(p_k^e)(S(\frac n{p_k^e},k)+[e\ne 1]) \]

\(\clubsuit\)最后：

将第一部分定义的积性函数 \(g(x)\) ，根据 \(g(p)=f(p)=p*(p-1)\) 的定义，拆成 \(g1(x)=x\) ，\(g2(x)=x^2\) 进行求解。

#include<bits/stdc++.h>
#define mem(a,b) memset(a,b,sizeof(a))
using namespace std;
typedef long long ll;
const int maxn=1e6+50;
const int mod=1e9+7;

const int inv6=166666668;
const int inv2=500000004;
int pri[maxn],ncnt,T,m,id1[maxn],id2[maxn];
bool isp[maxn];
ll n,a[maxn],g1[maxn],g2[maxn],sum1[maxn],sum2[maxn];
void init()
{
    T=sqrt(n+0.5);
    for(int i=2;i<=T;i++){
        if(!isp[i])pri[++ncnt]=i,sum1[ncnt]=(sum1[ncnt-1]+i)%mod,sum2[ncnt]=(sum2[ncnt-1]+1ll*i*i%mod)%mod;
        for(int j=1;j<=ncnt&&i*pri[j]<=T;j++){
            isp[i*pri[j]]=true;
            if(i%pri[j]==0)break;
        }
    }
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        a[++m]=n/l;
        if(a[m]<=T)id1[a[m]]=m;else id2[n/a[m]]=m;
        g1[m]=a[m]%mod;
        g2[m]=(g1[m]*(g1[m]+1)%mod*(g1[m]*2+1)%mod*inv6%mod-1+mod)%mod;
        g1[m]=(g1[m]*(g1[m]+1)%mod*inv2%mod-1+mod)%mod;
    }
    for(int i=1;i<=ncnt;i++){
        ll p2=1ll*pri[i]*pri[i],v;
        for(int j=1,id;j<=m&&p2<=a[j];j++){
            v=a[j]/pri[i];
            id=v<=T?id1[v]:id2[n/v];
            g1[j]=(g1[j]-1ll*pri[i]*(g1[id]-sum1[i-1]+mod)%mod+mod)%mod;
            g2[j]=(g2[j]-1ll*pri[i]*pri[i]%mod*(g2[id]-sum2[i-1]+mod)%mod+mod)%mod;
        }
    }
}
ll S(ll x,int y)
{
    if(pri[y]>=x)return 0;
    int id=x<=T?id1[x]:id2[n/x];
    ll sum=(g2[id]-g1[id]-(sum2[y]-sum1[y]+mod)%mod+mod+mod)%mod;
    for(int k=y+1;k<=ncnt&&1ll*pri[k]*pri[k]<=x;k++){
        ll pe=pri[k];
        for(int e=1;pe<=x;pe*=pri[k],e++)
            sum=(sum+pe%mod*((pe-1)%mod)%mod*(S(x/pe,k)+(e!=1))%mod)%mod;
    }
    return sum;
}
int main()
{
    scanf("%lld",&n);
    init();
    printf("%lld\n",(S(n,0)+1)%mod);
}

---

例题2：#6053.简单的函数

题意：

设 \(f(x)\) 有如下性质：

• \(f(1)=1\)
• \(f(p^c)=p\bigoplus c\) （\(p\) 是质数，\(\bigoplus\) 表示异或）
• \(f(ab)=f(a)f(b)\) （\(a\) 与 \(b\) 互质）

对于给定 \(n\) ，\(1\le n\le 10^{10}\) ，求 \(\sum^n_{i=1}f(i)\) 。对 \(10^9+7\) 取模。

解：

可以看出，除了 \(f(2)=3\) 外，对于 \(p\ne 2\) ，均有\(f(p)=p-1\) 。

同例题1，此时设 \(g1(x)=x\) ， \(g2(x)=1\) 。

则此计算到结果时，除了 \(f(1)\) 没有计算外，还将 \(f(2)\) 的结果计算成 \(1\) ，所以最后答案还得加上 \(3\) 。

#include<bits/stdc++.h>
#define mem(a,b) memset(a,b,sizeof(a))
using namespace std;
typedef long long ll;
const int maxn=1e6+50;
const int mod=1e9+7;

const int inv2=500000004;
int pri[maxn],ncnt,T,m,id1[maxn],id2[maxn];
bool isp[maxn];
ll n,a[maxn],g1[maxn],g2[maxn],sum1[maxn];
void init()
{
    T=sqrt(n+0.5);
    for(int i=2;i<=T;i++){
        if(!isp[i])pri[++ncnt]=i,sum1[ncnt]=(sum1[ncnt-1]+i)%mod;
        for(int j=1;j<=ncnt&&i*pri[j]<=T;j++){
            isp[i*pri[j]]=true;
            if(i%pri[j]==0)break;
        }
    }
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        a[++m]=n/l;
        if(a[m]<=T)id1[a[m]]=m;else id2[n/a[m]]=m;
        g1[m]=a[m]%mod;
        g1[m]=(g1[m]*(g1[m]+1)%mod*inv2%mod-1+mod)%mod;
        g2[m]=(a[m]-1+mod)%mod;
    }
    for(int i=1;i<=ncnt;i++){
        ll p2=1ll*pri[i]*pri[i],v;
        for(int j=1,id;j<=m&&p2<=a[j];j++){
            v=a[j]/pri[i];
            id=v<=T?id1[v]:id2[n/v];
            g1[j]=(g1[j]-1ll*pri[i]*(g1[id]-sum1[i-1]+mod)%mod+mod)%mod;
            g2[j]=(g2[j]-(g2[id]-(i-1)+mod)%mod+mod)%mod;
        }
    }
}
ll S(ll x,int y)
{
    if(pri[y]>=x)return 0;
    int id=x<=T?id1[x]:id2[n/x];
    ll sum=(g1[id]-g2[id]-(sum1[y]-y+mod)%mod+mod+mod)%mod;
    for(int k=y+1;k<=ncnt&&1ll*pri[k]*pri[k]<=x;k++){
        ll pe=pri[k];
        for(int e=1;pe<=x;pe*=pri[k],e++)
            sum=(sum+1ll*(pri[k]^e)%mod*(S(x/pe,k)+(e!=1))%mod)%mod;
    }
    return sum;
}
int main()
{
    scanf("%lld",&n);
    if(n==1){
        puts("1");return 0;
    }
    init();
    printf("%lld\n",(S(n,0)+3)%mod);
}

---

例题3：2020ccpc网络赛 Graph Theory Class

题意：

有 \(T\) 个测试用例，\(T\le 50\) 。 对于每个测试用例给定的 \(n,k\) ，输出 $\frac{n(n+3)}2+f(n+1)-4 $ 在模 \(k\) 下的结果。

其中 $f(x)=\sum_{1\le p\le x}p $。

解：

此时设 \(g(x)=x\) ，只需计算 \(G(n)\) 即为答案。

#include<bits/stdc++.h>
#define mem(a,b) memset(a,b,sizeof(a))
using namespace std;
typedef long long ll;
const int maxn=1e6+50;

int mod,inv2;
ll kpow(ll a,ll b){
    ll ans=1;
    while(b){
        if(b&1)ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
int pri[maxn],ncnt,T,m,id1[maxn],id2[maxn];
bool isp[maxn];
ll n,a[maxn],g1[maxn],sum1[maxn];
void init()
{
    m=ncnt=0;
    T=sqrt(n+0.5);
    for(int i=2;i<=T;i++){
        if(!isp[i])pri[++ncnt]=i,sum1[ncnt]=(sum1[ncnt-1]+i)%mod;
        for(int j=1;j<=ncnt&&i*pri[j]<=T;j++){
            isp[i*pri[j]]=true;
            if(i%pri[j]==0)break;
        }
    }
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        a[++m]=n/l;
        if(a[m]<=T)id1[a[m]]=m;else id2[n/a[m]]=m;
        g1[m]=a[m]%mod;
        g1[m]=(g1[m]*(g1[m]+1)%mod*inv2%mod-1+mod)%mod;
    }
    for(int i=1;i<=ncnt;i++){
        ll p2=1ll*pri[i]*pri[i],v;
        for(int j=1,id;j<=m&&p2<=a[j];j++){
            v=a[j]/pri[i];
            id=v<=T?id1[v]:id2[n/v];
            g1[j]=(g1[j]-1ll*pri[i]*(g1[id]-sum1[i-1]+mod)%mod+mod)%mod;
        }
    }
}

int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%d",&n,&mod);
        inv2=kpow(2,mod-2);
        n++;init();n--;n%=mod;
        ll x=g1[id2[1]];
        printf("%lld\n",(x-4+mod+n%mod*(n+3)%mod*inv2%mod)%mod);
    }
}