LOJ

分析

注意到

$$ f(p)=\begin{cases}p-1,&p\neq 2\\p+1,&p=2\end{cases} $$

不妨假装 $f(2)=1$,于是就可以 Min_25 筛了。

考虑到在 Min_25 筛的第二步即求 $S(n,j)$ 时,如果 $j=1$ 则质数部分算了 $2$ 的贡献,直接加上 $2$ 即可。

代码

// ===================================
//   author: M_sea
//   website: http://m-sea-blog.com/
// ===================================
#include <bits/stdc++.h>
#define re register
using namespace std;
typedef long long ll;

inline ll read() {
    ll X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X*w;
}

const int N=200000+10;
const int mod=1e9+7,inv2=500000004;

ll n; int s;

int p[N],v[N],cnt=0;
inline void sieve(int n) {
    for (re int i=2;i<=n;++i) {
        if (!v[i]) p[++cnt]=i;
        for (re int j=1;j<=cnt&&i*p[j]<=n;++j) {
            v[i*p[j]]=1; if (i%p[j]==0) break;
        }
    }
}

ll w[N]; int id1[N],id2[N],top=0;
int g0[N],g1[N],gs1[N];

inline int S(ll x,int j) {
    if (x<=1||p[j]>x) return 0;
    int k=x<=s?id1[x]:id2[n/x];
    int res=((1ll*g1[k]-g0[k]-gs1[j-1]+(j-1))%mod+mod)%mod;
    if (j==1) res+=2;
    for (re int i=j;i<=cnt&&1ll*p[i]*p[i]<=x;++i) {
        ll p1=p[i],p2=1ll*p[i]*p[i];
        for (re int e=1;p2<=x;++e,p1=p2,p2*=p[i])
            res=(res+1ll*S(x/p1,i+1)*(p[i]^e)%mod+(p[i]^(e+1))%mod)%mod;
    }
    return res;
}

int main() {
    n=read(),s=sqrt(n); sieve(s);
    for (re int i=1;i<=cnt;++i) gs1[i]=(gs1[i-1]+p[i])%mod;
    for (re ll l=1,r;l<=n;l=r+1) {
        r=n/(n/l); w[++top]=n/l;
        n/l<=s?id1[n/l]=top:id2[n/(n/l)]=top;
        int x=(n/l)%mod;
        g0[top]=(x-1+mod)%mod;
        g1[top]=(1ll*x*(x+1)%mod*inv2%mod-1+mod)%mod;
    }
    for (re int i=1;i<=cnt;++i)
        for (re int j=1;j<=top&&1ll*p[i]*p[i]<=w[j];++j) {
            ll x=w[j]/p[i]; int k=x<=s?id1[x]:id2[n/x];
            g0[j]=(g0[j]-(g0[k]-(i-1)+mod)%mod+mod)%mod;
            g1[j]=(g1[j]-1ll*p[i]*(g1[k]-gs1[i-1]+mod)%mod+mod)%mod;
        }
    printf("%d\n",(S(n,1)+1)%mod);
    return 0;
}
最后修改:2020 年 02 月 25 日 01 : 43 PM