LOJ 6053 简单的函数 - Min_25筛

不知道为啥脑子一抽就来温习了一波Min_25筛
这题F( p )=p-1（除了p=2，这个特判一下即可），最后把2加回来。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define Rep(i,v) rep(i,0,(int)v.size()-1)
#define lint long long
#define mod 1000000007
#define ull unsigned lint
#define db long double
#define pb push_back
#define mp make_pair
#define fir first
#define sec second
#define gc getchar()
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
typedef pair<int,int> pii;
typedef set<int>::iterator sit;
#define P(x) (x>=mod?x-=mod:0)
#define Q(x) (x<0?x+=mod:0)
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%mod) (k&1)?ans=(lint)ans*x%mod:0;return ans; }
namespace min_25{
	const int MXS=100010;lint n;
	int s,A[MXS],B[MXS],pri[MXS],cnt;
	int a[MXS],b[MXS],mip[MXS],np[MXS];
	inline int MySqrt(lint n)
	{
		int x=(int)sqrt(n);
		while((lint)x*x<n) x++;
		while((lint)x*x>n) x--;
		return x;
	}
	inline int init()
	{
		memset(np,0,sizeof(int)*(s+1));
		memset(A,0,sizeof(int)*(s+1));
		memset(B,0,sizeof(int)*(s+1));
		return cnt=0;
	}
	inline int prelude()
	{
		rep(i,2,s) if(!np[i]) rep(j,i,s/i) np[i*j]=1;
		rep(i,2,s) if(!np[i]) pri[++cnt]=i;return pri[++cnt]=s+1;
	}
	inline int getg(int xs)
	{
	//	rep(i,1,s) debug(i)sp,debug(a[i])sp,debug(b[i])ln;cerr ln;
		rep(p,2,s) if(!np[p])
		{
			lint r=(lint)p*p;int ed1=s/p,ed2=(int)min(n/r,(lint)s);
			rep(i,1,ed1) b[i]-=(lint)mip[p]*(b[i*p]-a[p-1])%mod,P(b[i]),Q(b[i]);
			rep(i,ed1+1,ed2) b[i]-=(lint)mip[p]*(a[n/i/p]-a[p-1])%mod,P(b[i]),Q(b[i]);
			for(int i=s;i>=r;i--) a[i]-=(lint)mip[p]*(a[i/p]-a[p-1])%mod,P(a[i]),Q(a[i]);
	//		debug(p)ln;rep(i,1,s) debug(i)sp,debug(a[i])sp,debug(b[i])ln;cerr ln;
		}
		rep(i,1,s) A[i]=(A[i]+(lint)a[i]*xs)%mod,B[i]=(B[i]+(lint)b[i]*xs)%mod;return 0;
	}
	inline int F(int p,int c) { return p^c; }
	inline int G(lint x) { return x<=s?A[x]:B[n/x]; }
	int S(lint x,int y,int curd=0)
	{
	//	debug(x)sp,debug(y)sp,debug(pri[y])ln;
		if(pri[y]>x) return 0;
		int res=G(x)-G(pri[y]-1);Q(res);
		rep(i,y,cnt&&(lint)pri[i]*pri[i]<=x)
		{
			lint v=pri[i];
			for(int j=1;v<=x/pri[i];v*=pri[i],j++)
				res=(res+(lint)F(pri[i],j)*S(x/v,i+1,curd+1)+F(pri[i],j+1))%mod;
		//		debug(pri[i])sp,debug(j)sp,debug(j+1)ln,
		//		debug(qwq)sp,debug(F(pri[i],j))sp,debug(F(pri[i],j+1))ln ln;
		}
	//	debug(x)sp,debug(pri[y])sp,debug(res)ln ln;
		return res;
	}
	inline int solve(lint _n)
	{
		s=MySqrt(n=_n),init(),prelude();
		rep(i,2,s) if(!np[i]) mip[i]=1;
		rep(i,1,s) a[i]=i,b[i]=(n/i)%mod;getg(mod-1);
	//	rep(i,1,s) debug(i)sp,debug(A[i])sp,debug(B[i])ln;
		
		rep(i,2,s) if(!np[i]) mip[i]=i;int inv2=fast_pow(2,mod-2),x;
		rep(i,1,s) a[i]=(i*(i+1ll)/2)%mod,x=(n/i)%mod,b[i]=x*(x+1ll)%mod*inv2%mod;
		getg(1);
	//	rep(i,1,s) debug(i)sp,debug(A[i])sp,debug(B[i])ln;cerr ln;
		int Ans=S(n,1)+1+(n>=2?2:0);return P(Ans),Q(Ans),Ans;
	}
}
int main()
{
	lint n;cin>>n;return printf("%d\n",min_25::solve(n));
}