【洛谷P5325】【模板】Min_25筛

最新推荐文章于 2022-01-20 22:21:31 发布

原创最新推荐文章于 2022-01-20 22:21:31 发布 · 577 阅读

1 ·

CC 4.0 BY-SA版权

积性函数同时被 2 个专栏收录

4 篇文章

订阅专栏

min_25筛

2 篇文章

订阅专栏

博客提及n的范围小于1010，建议抛弃ees和洲阁筛，选择min_25，并给出了相关代码，主要围绕算法选择展开。

在这里插入图片描述
n的范围是小于 $10^{10}$

抛弃ees和洲阁筛，拥抱min_25
代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<vector>
#include<algorithm>
#define INF 0x7f7f7f7f
#define maxx 100050
#define mod 1000000007
using namespace std;
typedef long long ll;
ll n,N;
int pr[maxx],cnt;
bool isp[maxx];
int p1[maxx],p2[maxx];
void init()
{
    N=sqrt(n);
    while((N+1)*(N+1)<=n)N++;
    for(int i=2;i<=N;i++)
    {
        if(!isp[i])
        {
            pr[++cnt]=i;
            p1[cnt]=(p1[cnt-1]+i)%mod;
            p2[cnt]=(p2[cnt-1]+(ll)i*i)%mod;
        }
        for(int j=1;j<=cnt&&(ll)i*pr[j]<=N;j++)
        {
            isp[i*pr[j]]=true;
            if(i%pr[j]==0)break;
        }
    }
}
ll w[maxx<<1],tot;
int g1[maxx<<1],g2[maxx<<1];
int id1[maxx],id2[maxx];
int inv2=500000004,inv6=166666668;
void G()
{
    for(ll i=1,last;i<=n;i=last+1)
    {
        ll now=n/i;
        w[++tot]=now;
        last=n/now;
        if(now<=N)id1[now]=tot;
        else id2[last]=tot;
    }
    for(int i=1;i<=tot;i++)
    {
        ll temp=w[i]%mod;
        g1[i]=temp*(temp+1)%mod*inv2%mod-1;
        g2[i]=(temp+1)*temp%mod*(2*temp+1)%mod*inv6%mod-1;
    }

    for(int i=1;i<=cnt;i++)
    {
        ll _p=(ll)pr[i]*pr[i]%mod;
        for(int j=1;j<=tot&&(ll)pr[i]*pr[i]<=w[j];j++)
        {
            ll now=w[j]/pr[i];
            int k=now<=N?id1[now]:id2[n/now];
            g1[j]-=((ll)g1[k]-p1[i-1]+mod)*pr[i]%mod;
            if(g1[j]<0)g1[j]+=mod;

            g2[j]-=(g2[k]-p2[i-1]+mod)*_p%mod;
            if(g2[j]<0)g2[j]+=mod;
        }
    }
}
int S(ll x,int y)
{
    if(x<=pr[y])return 0;
    int k=x<=N?id1[x]:id2[n/x];
    int ans=(g2[k]-g1[k]+mod)%mod;
    ans-=(p2[y]-p1[y]+mod)%mod;
    if(ans<0)ans+=mod;

    for(int i=y+1;i<=cnt&&(ll)pr[i]*pr[i]<=x;i++)
    {
        ll pe=pr[i];
        for(int e=1;pe<=x;e++,pe*=pr[i])
        {
            ll _pe=pe%mod;
            ans+=_pe*(_pe-1+mod)%mod*(S(x/pe,i)+(e!=1))%mod;
            ans%=mod;
        }
    }
    return ans;
}
int main()
{
    cin>>n;
    init();
    G();
    //cout<<g1[1]<<endl;
    cout<<(S(n,0)+1)%mod<<endl;
    return 0;
}