【LOJ6053】—简单的函数（min_25筛)

传送门

考虑对于一个质数
$f(p)= egin{cases} p-1 && p>2 \ p+1 && p=2\ end{cases}$

所以求出$F(p)=1$和$s(p)=p$
相减就可以了求出$f$了

剩下的就一样了

不会$min_{25}$的看这个

#include<bits/stdc++.h>
using namespace std;
#define gc getchar
inline int read(){
	char ch=gc();
	int res=0,f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
	return f?res:-res;
}
#define re register
#define pb push_back
#define cs const
#define pii pair<int,int>
#define fi first
#define se second
#define ll long long
const int N=1000006;
const int mod=1e9+7;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline void Dec(int &a,int b){(a-=b)<0?(a+=mod):0;}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;
}
const int inv2=ksm(2,mod-2);
int f1[N],f2[N],s1[N],s2[N],pr[N];
int lim,tot;
inline int C(ll x){
	return mul(mul((x%mod),(x+1)%mod),inv2);
}
inline void init(ll k){
	if(k<=1)return ;
	lim=sqrt(k);tot=0;
	for(int i=1;i<=lim;i++)f1[i]=i-1,f2[i]=(k/i-1)%mod,s1[i]=dec(C(i),1),s2[i]=dec(C(k/i),1);
	for(int p=2;p<=lim;p++){
		if(f1[p]==f1[p-1])continue;
		pr[++tot]=p;
		for(int i=1;i<=lim/p;i++)Dec(f2[i],dec(f2[i*p],f1[p-1])),Dec(s2[i],mul(p,dec(s2[i*p],s1[p-1])));
		for(int i=lim/p+1;1ll*i*p*p<=k&&i<=lim;i++)
			Dec(f2[i],dec(f1[k/i/p],f1[p-1])),Dec(s2[i],mul(p,dec(s1[k/i/p],s1[p-1])));
		for(int i=lim;i>=1ll*p*p;i--)Dec(f1[i],dec(f1[i/p],f1[p-1])),Dec(s1[i],mul(p,dec(s1[i/p],s1[p-1])));
	}
	for(int i=1;i<=lim;i++){
		Dec(s1[i],f1[i]),Add(s1[i],mul(2,i>=2));
		Dec(s2[i],f2[i]),Add(s2[i],mul(2,(k/i)>=2ll));
	}
}
ll n;
int ans;
void calc(int pos,int coef,ll res){
	int g=(res<=lim?s1[res]:s2[n/res]);
	Add(ans,mul(coef,dec(g,s1[pr[pos-1]])));
	for(int i=pos;i<=tot;i++){
		if(1ll*pr[i]*pr[i]>res)return;
		for(ll now=pr[i],xs=1;now<=res;now=now*pr[i],xs++){
			if(now*pr[i]<=res)calc(i+1,mul(coef,pr[i]^xs),res/now);
			if(xs>=2)Add(ans,mul(coef,pr[i]^xs));
		}
	}
}
signed main(){
	scanf("%lld",&n);
	ans=1,init(n);
	calc(1,1,n);
	cout<<ans<<'
';
}