前置知识:min25 筛,即你要用 min25 筛通过板子题,不管写成啥样,不管常数多大,但是你要了解一点 min25。
没有特殊说明的话有如下记号(大部分记号与 oi-wiki 一致):
- \(x/y=\lfloor\frac{x}{y}\rfloor\)
- \(\text{P}\) 为素数集合,\(p_k\) 表示第 \(k\) 小素数。特别地,令 \(p_{0} = 1\)。\(\text{lpf}(n)\) 表示 \(n\) 的最小素因子。
min25 筛说的是这样一个东西,记:
\(F_{\mathrm{pr}}(n) := \sum\limits_{2 \le p \le n,p\in \text{P}} f(p)\)
\(F_{k}(n) := \sum\limits_{2\le i\le n} [p_{k} \le \operatorname{lpf}(i)] f(i)\)
则 \(F_{k}(n)= \sum\limits_{\substack{k \le i \\ p_{i}^{2} \le n}} \sum\limits_{\substack{c \ge 1 \\ p_{i}^{c + 1} \le n}} \left(f\left(p_{i}^{c}\right) F_{i + 1}\left(n / p_{i}^{c}\right) + f\left(p_{i}^{c + 1}\right)\right) + F_{\mathrm{pr}}(n) - F_{\mathrm{pr}}(p_{k - 1})\)
接下来先上常数小的代码:
$\texttt{code}$
//落谷 P5325
//https://www.luogu.com.cn/problem/P5325
#include<bits/stdc++.h>
#define u32 unsigned int
#define u64 unsigned long long
#define u128 unsigned __int128
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const u32 N=2e5+5,M=227649,mod=1e9+7,inv2=(mod+1)/2,inv6=(mod+1)/6;
u32 m,cnt,pr[N],g[N][2],tot,vp[M][45];
u64 n,to[N],inv[N];
inline u64 div(u64 x,u32 y){return (u128)x*inv[y]>>64;}
inline constexpr u32 AD(u32 x,u32 y){return x+=y,x>=mod?x-mod:x;}
inline constexpr u32 S1(u64 x){x%=mod;return (u64)x*(x+1)%mod*inv2%mod;}
inline constexpr u32 S2(u64 x){x%=mod;return (u64)x*(x+1)%mod*(x<<1|1)%mod*inv6%mod;}
inline constexpr u32 ff(u64 x){x%=mod;return x*(x-1)%mod;}
inline u32 TO(u64 x){return x<=m?x:tot-(u64)((double)(n)/x)+1;}
u32 S(u64 x,u32 y)
{
if(pr[y]>=x) return 0;u32 ans=g[TO(x)][0]+mod-g[pr[y]][0];u64 X;
for(u32 i=y+1,e;i<=cnt&&(u64)pr[i]*pr[i]<=x;i++)
for(e=1,X=div(x,pr[i]);X>=pr[i];e++,X=div(X,pr[i]))
ans=(ans+vp[i][e+1]+(u64)vp[i][e]*S(X,i))%mod;
return ans;
}
int main()
{
scanf("%lld",&n);m=sqrt(n);
for(u32 i=1;i<=m;i++) inv[i]=~0ull/i+1;
for(u64 i=1,t=n;i<=n;i=n/t+1,t=n/i) to[++tot]=n/t,g[tot][0]=S2(to[tot])-1,g[tot][1]=S1(to[tot])-1;//一定要减去1的贡献!
for(u32 i=2;i<=m;i++) if(g[i][0]^g[i-1][0])
{
pr[++cnt]=i;u32 t0=g[i-1][0]+mod,t1=g[i-1][1]+mod,pf=(u64)i*i%mod;
for(u32 j=tot,t=TO(div(to[j],i));t>=i;t=TO(div(to[--j],i)))
g[j][0]=(g[j][0]+(u64)pf*(t0-g[t][0]))%mod,
g[j][1]=(g[j][1]+(u64)i*(t1-g[t][1]))%mod;
}
for(u32 i=1;i<=tot;i++) g[i][0]=AD(g[i][0],mod-g[i][1]);
for(u32 i=1;i<=cnt;i++) for(u64 s=pr[i],e=1;s<=n*pr[i];s*=pr[i],e++) vp[i][e]=ff(s);
printf("%d",AD(S(n,0),1));
return 0;
}//n<=10^12 内不丢精度