P5325 【模板】Min_25筛(min25筛)

P5325 【模板】Min_25筛(min25筛)

题目大意

定义积性函数f(x),且 f ( p k ) = p k ( p k − 1 ) f(p^k)=p^k(p^k-1) f(pk)=pk(pk1)(pp是一个质数),求

∑ i = 1 n f ( x ) ​ \sum_{i=1}^n f(x)​ i=1nf(x)

对10^9+7取模。

解题思路

min25筛版题https://blog.csdn.net/baiyifeifei/article/details/90454317

AC代码

// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;
#define int long long 
const int mod=1e9+7;
const int size=2e6+5;
typedef long long LL;
int tot=0,tol=0;
LL num[size];
bool prime[size];int  p[size];
LL h1[size],h2[size];
LL pre1[size],pre2[size];
LL n; 
LL s;
inline int ID(LL nums){if(nums<=s) return nums;return tot-n/nums+1;}
int inv2,inv3,inv6;
LL quick_pow(LL a,LL b)
{
    LL ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return ans;
}
LL get_g(LL n,int k)
{
    if(n<=1||p[k]>n) return 0;
    int idn=ID(n);
    LL ans=(h2[idn]-h1[idn]+mod)%mod;
    ans=(ans-(pre2[k-1]-pre1[k-1])+mod)%mod;
    while(1LL*p[k]*p[k]<=n&&k<=tol)
    {
        LL pk=p[k],t=p[k];
        for(int e=1;t<=n;e++,t=t*pk)
        {
            ans=(ans+(t%mod)*(t%mod-1)%mod*((e>1)+get_g(n/t,k+1))%mod)%mod;
        }
        k++;
    }
    return ans%mod;
}
void get_h()
{
    s=sqrt(n);
    pre1[0]=pre2[0]=0;
    while(1LL*s*s<=n) s++;
    while(1LL*s*s>n) s--;
    for(int i=1;i<=s;i++) prime[i]=true;
    for(int i=2;i<=s;i++)
    {
        if(prime[i])
        {
            p[++tol]=i;
            pre1[tol]=(pre1[tol-1]+i)%mod;
            pre2[tol]=(pre2[tol-1]+1LL*i*i%mod)%mod;
        }
        for(int j=1;j<=tol&&p[j]*i<=s;j++)
        {
            prime[p[j]*i]=false;
            if(i%p[j]==0) break;
        } 
    }
    tot=0;
    for(int i=1;i<=s;i++) num[++tot]=i;
    for(int i=s;i>=1;i--) if(n/i>s) num[++tot]=n/i;
    for(int i=1;i<=tot;i++) h1[i]=h2[i]=num[i]%mod;//num最大可以达到n需要做这一步,否则会溢出LL 
    for(int i=1;i<=tot;i++) h2[i]=1LL*h2[i]*(h2[i]+1)%mod*(2*h2[i]+1)%mod*inv6%mod-1;//计算的时候将1减去 
    for(int i=1;i<=tot;i++) h1[i]=1LL*(1+h1[i])*(h1[i])%mod*inv2%mod-1;
    int x=1;
    for(int j=1;j<=tol;j++)
    {
        while(num[x]<1LL*p[j]*p[j]) x++;
        for(int i=tot;i>=x;i--) 
        {
            h1[i]=((h1[i]-1LL*p[j]*((h1[ID(num[i]/p[j])]-pre1[j-1]+mod)%mod)%mod)%mod+mod)%mod;
            h2[i]=((h2[i]-(1LL*p[j]*p[j]%mod*(h2[ID(num[i]/p[j])]-pre2[j-1]+mod)%mod)%mod)%mod+mod)%mod;
        }
    }
}		
int32_t main()
{
    inv2=quick_pow(2,mod-2),inv3=quick_pow(3,mod-2);
    inv6=1LL*inv2*inv3%mod;
    scanf("%lld",&n);
    get_h();
    printf("%lld\n",(get_g(n,1)+1)%mod);
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值