Min25 sieve

Defining a product function $f(x)$,$f(p^c)$ It can be expressed as multiple completely integrable functions $h_i(p^c)$ Add up , seek
$$\sum _{i=1}^{n}f(i)$$

$$\sum _{i=1}^{n}f(i)=\sum _{p\le n}f(p)+\sum _{iisnotaprime}f(i)$$
set up $g(i,j,k)$ Express ${\sum }_{w=1}^n[w\in prime|minp>p_j]h_k(w)$
$$g(i,j,k)=g(i,j-1,k)-[g(\lfloor \frac{i}{p_j}\rfloor ,j-1,k)-g(p_j-1,j-1,k)]h_k(p_j)$$
$g(i,maxprime,k)$ said ${\sum }_{p\le i}{h}_k(p)$

set up $S(n,x)$ Express ${\sum }_{w=1}^n[minp>p_j]f(w)$
$$S(n,x)=H(n)-H(p_x)+\sum _{p_k^e\le n,k>x}[S(\lfloor \frac{n}{p^e}\rfloor ,k)+(e\ne 1)]f(p^e)$$
${\sum }_{i=1}^{n}f(i)=S(n,0)$
Fine implementation , It can be done $O(\frac{n^{\frac{3}{4}}}{\log n})$
P5325 【 Templates 】Min_25 sieve

#include<cstdio>
#include<cmath>
#define ll long long
using namespace std;
const ll mod=1e9+7,inv2=500000004,inv3=333333336;
ll ksm(ll a,ll b)
{
    
	ll res=1;
	while(b)
	{
    
		if(b&1) res=res*a%mod;
		a=a*a%mod,b>>=1;
	}
	return res;
}
const int N=2e6+1000;
ll n,w[N+10],sp1[N+10],sp2[N+10],g1[N+10],g2[N+10];int m,cnt,id1[N+10],id2[N+10],prime[N+10];
bool vis[N+10];
int id(ll x)
{
    
	if(x<=m) return id1[x];
	else return id2[n/x];
}
void init()
{
    
	m=sqrt(n);
	for(int i=2;i<=m;i++)
	{
    
		if(!vis[i]) prime[++prime[0]]=i,sp1[prime[0]]=(sp1[prime[0]-1]+i)%mod,sp2[prime[0]]=(sp2[prime[0]-1]+1ll*i*i%mod)%mod;
		for(int j=1;j<=prime[0]&&i*prime[j]<=m;j++)
		{
    
			vis[i*prime[j]]=true;
			if(i%prime[j]==0) break;
		}
	}
}
ll S(ll n,int x)
{
    
	if(prime[x]>=n) return 0;
	int l=id(n);
	ll res=(g2[l]+mod-g1[l]+mod-(sp2[x]-sp1[x]))%mod;
	for(int i=x+1;i<=prime[0]&&1ll*prime[i]*prime[i]<=n;i++)
	{
    
		ll pe=prime[i],tmp;
		for(int e=1;pe<=n;e++,pe=pe*prime[i]) tmp=pe%mod,res=(res+tmp*(tmp-1)%mod*(S(n/pe,i)+(e!=1))%mod)%mod; 
	}
	return res;
}
int main()
{
    
	scanf("%lld",&n),init();
	for(ll l=1,r;l<=n;l=r+1)
	{
    
		r=n/(n/l),w[++cnt]=n/l;
		g1[cnt]=(w[cnt]%mod+2)*(w[cnt]%mod+mod-1)%mod*inv2%mod,g2[cnt]=(w[cnt]%mod*(w[cnt]%mod+1)%mod*(w[cnt]*2%mod+1)%mod*inv2%mod*inv3%mod+mod-1)%mod;
		if(w[cnt]<=m) id1[w[cnt]]=cnt;
		else id2[n/w[cnt]]=cnt;
	}
	for(int i=1;i<=prime[0];i++)
		for(int j=1;j<=cnt&&1ll*prime[i]*prime[i]<=w[j];j++)
			g1[j]=(g1[j]+mod-(g1[id(w[j]/prime[i])]+mod-sp1[i-1])%mod*prime[i]%mod)%mod,
			g2[j]=(g2[j]+mod-(g2[id(w[j]/prime[i])]+mod-sp2[i-1])%mod*prime[i]%mod*prime[i]%mod)%mod;
	printf("%lld",(S(n,0)+1)%mod);
	return 0;
}