LOJ 6053 簡單的函式 - Min_25篩
阿新 • • 發佈:2018-12-30
不知道為啥腦子一抽就來溫習了一波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));
}