一个筛子(例一 LOJ 6053)

一个筛子(例一 LOJ 6053)

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cstring>
#include <vector>
#include <map>
#include <queue>
#include <iterator>
#include <utility>
#include <algorithm>
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;
typedef pair<int, int> pii;
#define mp make_pair
#define fi first
#define se second
#define All(x) (x).begin(), (x).end()
#define Y1 "YES"
#define N1 "NO"
#define ENDL '\n'
#define count2(x) __builtin_popcount(x)
#define countleadingzero(x) __builtin_clz(x)
inline ll read() {  // not solve LLONG_MIN LMAX=9,223,372,036,854,775,807
    ll s = 0, w = 1;
    char ch = getchar();
    while (ch < '0' || ch > '9') {
        if (ch == '-')
            w = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9') s = s * 10 + ch - '0', ch = getchar();
    return s * w;
}
struct ygsz{
	typedef __int128_t returnT;
	typedef long long argT;
	typedef unsigned int sqargT;
	static constexpr argT maxm=1e5+10;
	sqargT primelist[maxm],pcnt=0;
	bool isnprime[maxm];
	double inv[maxm];
	static constexpr argT mod = 1e9+7;
	inline argT fpowll(argT a,argT b,argT p){//a^b mod p;
		assert(p!=0);
		argT ans=1%p,base=a%p;
		for(;b>0;b/=2){
			if(b&1)ans=(returnT)ans*base%p;
			base=(returnT)base*base%p;
		}
		return ans;
	}
	const argT inv2=fpowll(2,mod-2,mod);
	inline void init(){
	    constexpr argT m=maxm-1;
	//    constexpr argT mz=m/log10(m);
	    primelist[++pcnt]=2;
	    for(argT i = 1;i <= m; ++i)
	        inv[i] = 1.0 / i;
	    // i*i may overflow
	    for(argT j = 4; j <=m; j+=2)
	    isnprime[j] = true;
	    for (argT i = 3; i <= m; i+=2){
	        if(!isnprime[i]){
	            primelist[++pcnt]=i;
	            // i*i may overflow
	            for(argT j = i * i; j <=m; j+=i)
	            isnprime[j] = true;
	        }
	    }
	    return ;
	}
	argT a[maxm];
	returnT b[maxm];//质数前缀和
	argT s1[maxm],s2[maxm];//质数个数前缀和 
	inline argT f(const argT x){
		return x;
	}
	//int ppt=0;
	inline argT n34divlognsieve(const argT n) {
	    if (n <= 1) return 0;
	    const sqargT m = sqrt(n);
	    constexpr double eps = 1e-7;
	    for (sqargT i = 2; i <= m; ++i)    
	        a[i] = (a[i - 1] + i)%mod,s1[i] = i - 1;
	    for (sqargT i = 1; i <= m; ++i) {  
	        returnT sum_n = n / i;
	        b[i] = ((sum_n + sum_n * sum_n) *inv2 - 1)%mod;
	        s2[i] = sum_n - 1;
	    }
	    const sqargT pt=upper_bound(primelist+1,primelist+1+pcnt,m)-primelist;
	    for (sqargT r=1;r<pt;++r) {
	        const sqargT prime=primelist[r];
	        const argT t = a[prime - 1];
	        const argT s = s1[prime - 1];
	        const double invprime=inv[prime];
	        const sqargT limmp = m * invprime + eps;
	        for (sqargT i = 1; i <= limmp; ++i){
	        	b[i] =(b[i]- (b[i * prime] - t) * prime)%mod; 
	        	s2[i] =(s2[i] - (s2[i *prime] - s))%mod;
				s2[i] =(s2[i] + mod)%mod;
			}
	            
	        const sqargT limnpp = n / prime / prime;
	        const argT limnp = n / prime;
	        const sqargT lim2=min(limnpp,m);
	        double *ptr=&inv[limmp + 1];
	        for (sqargT i = limmp + 1; i <= lim2 ; ++i,++ptr) {
	        	b[i] =(b[i]- (a[(unsigned int)(limnp * (*ptr)  + eps)] - t) * prime)%mod;
	        	s2[i] =(s2[i]- (s1[(unsigned int)(limnp * (*ptr)  + eps)] - s) )%mod;
	        	s2[i] =(s2[i]+mod)%mod;
			}
	            
	        const argT limpp = prime * prime;
	        for (argT i = m; i >= limpp; --i){
	        	a[i] =(a[i]- (a[(unsigned int)(i * invprime + eps)] - t) * prime)%mod;
	        	s1[i] =(s1[i] - (s1[(unsigned int)(i * invprime + eps)] - s))%mod;
	        	s1[i] =(s1[i]+mod)%mod;
			}
	    }
	    for (argT i = 1; i <= m; ++i){
	    	a[i]=(a[i]-s1[i]+(i>=2)*2+mod)%mod;
	    	b[i]=(b[i]-s2[i]+(n/i>=2)*2+mod)%mod;
		}
	    return b[1];
	}
	argT n;
	returnT work(argT pos,argT res,argT rem){
		argT t = (rem <= sqrt(n)? a[rem]: b[n/rem]);
		returnT ans = (res * (t -a[primelist[pos-1]])%mod)%mod;
		for(argT i = pos; i <= pcnt; ++i){
			if(primelist[i]*primelist[i]>rem)return ans;
			for(argT t=primelist[i],js=1;t<=rem;t*=primelist[i],++js){
				if(t*primelist[i]<rem)ans+=work(i+1,res*(primelist[i]^js)%mod,rem/t);
				if(js>=2)ans=(ans+res*(primelist[i]^js)%mod )%mod;
			}
		}
		return ans;
	}	
}P;

int main() {
    P.init();
    ll n=read();
    P.n=n;
   	P.n34divlognsieve(n);
    cout<<(int)P.work(1,1,n)+1;
    return 0;
}