Loj 6053(EES筛法)

problem

思路

满足EES筛法的要求
考虑 f ( p ) f(p) f(p) p 1 p^1 p1是啥,当 p > 2 p>2 p>2时, f ( p ) = p − 1 f(p)=p-1 f(p)=p1,而 f ( 2 ) = 3 f(2)=3 f(2)=3
所以预处理出 p 1 , p 0 p^1,p^0 p1,p0的前缀和即可
并注意,当 i > = 2 i>=2 i>=2时,这个前缀和少计算了2(即p=2的贡献3-1),所以加上一下

在dfs中维护时,没什么特别的,先计数第一类贡献(当前数乘以更大的素数的贡献),再计数第二类贡献(最大的质因子有更大的指数)

注意别忘了取模

//phi(1^2)=1
//phi(p^2)=p^2-p
//phi((p^e)^2)=phi((p^{e-1})^2)*p*p;
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
const int ni6=166666668;
const int ni2=500000004;
ll n,M;
vector<int> pre[2],hou[2],primes;

inline int add(const int x, const int v) {
    return x + v >= mod ? x + v - mod : x + v;
}
inline int dec(const int x, const int v) {
    return x - v < 0 ? x - v + mod : x - v;
}

//这里res是n/枚举的数
int dfs(ll res, int last, ll f){
    //最大质因子是prime[last-1] 但将1放在外面值显然一样
    int t=dec((res > M ? hou[1][n/res] : pre[1][res]),pre[1][primes[last]-1]);
    int ret= t*f%mod;//这里需修改
    for(int i=last;i<(int) primes.size();++i){
        int p = primes[i];
        if((ll)p*p > res) break;
        int e=1;
        for(ll q=p,nres=res,nf=f*(p^e)%mod;q*p<=res;q*=p){//nf需修改
            ret = add(ret,dfs(nres/=p,i+1,nf));//枚举更大的数
            e++;
            nf = f* (p^e)%mod;
            ret =add(ret,nf);//指数大于1时,记上贡献
        }
    }
    return ret;
}

inline int ff(ll x){
    x%=mod;
    return x*(x+1)%mod*ni2%mod;
}

int solve(ll n){
    M=sqrt(n);
    for(int i=0;i<2;++i){
        pre[i].clear();pre[i].resize(M+1);
        hou[i].clear();hou[i].resize(M+1);
    }
    primes.clear();primes.reserve(M+1);
    for(int i=1;i<=M;++i){
        pre[0][i]=i-1;
        hou[0][i]=(n/i-1)%mod;
        pre[1][i]=dec(ff(i),1);
        hou[1][i]=dec(ff(n/i),1);
    }
    for(int p=2;p<=M;++p){
        if(pre[0][p]==pre[0][p-1]) continue;
        primes.push_back(p);
        const ll q=(ll)p*p,m=n/p;
        const int pnt0=pre[0][p-1],pnt1=pre[1][p-1];
        const int mid=M/p;
        const int End=min((ll)M,n/q);
        for(int i=1;i<=mid;++i){
            hou[0][i]=dec(hou[0][i],dec(hou[0][i*p],pnt0));
            hou[1][i]=dec(hou[1][i],dec(hou[1][i*p],pnt1)*(ll)p%mod);
            //hou[2][i]=dec(hou[2][i],dec(hou[2][i*p],pnt2)*q%mod);
        }
        for(int i=mid+1;i<=End;++i){
            hou[0][i]=dec(hou[0][i],dec(pre[0][m/i],pnt0));
            hou[1][i]=dec(hou[1][i],dec(pre[1][m/i],pnt1)*(ll)p%mod);
            //hou[2][i]=dec(hou[2][i],dec(pre[2][m/i],pnt2)*q%mod);
        }
        for(int i=M;i>=q;--i){
            pre[0][i]=dec(pre[0][i],dec(pre[0][i/p],pnt0));
            pre[1][i]=dec(pre[1][i],dec(pre[1][i/p],pnt1)*(ll)p%mod);
            //pre[2][i]=dec(pre[2][i],dec(pre[2][i/p],pnt2)*q%mod);
        }
    }
    primes.push_back(M+1);
    for (int i = 1; i <= M;++i) {
        pre[1][i] = i>=2 ? add(2,dec(pre[1][i], pre[0][i])) : (dec(pre[1][i], pre[0][i]));//p-1
        hou[1][i] = n>=2*i ? add(2,dec(hou[1][i], hou[0][i])) :(dec(hou[1][i], hou[0][i]));
    }
    return n>1 ? add(dfs(n,0,1),1) : 1;//别忘了加上1的贡献
}

int main()
{
    //freopen("in.txt","r",stdin);
    ios::sync_with_stdio(false);
    cin>>n;
    cout<<solve(n)<<endl;
    return 0;
}

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值