Min_25 筛用来求一类特殊数论函数的前缀和。
要求的函数 (f) 必须满足以下条件：
- (f) 是积性函数。
- 对于质数 (p)，(f(p)) 是一个关于 (p) 的简单低次多项式。
- 对于质数 (p) 和正整数 (e)，(f(p^e)) 可以快速计算。
转化题目
题目中的函数 (f(p^k)=p^k(p^k-1)) 显然满足要求。
即对于 (x=p^k)，(f(x)=x^2-x)，其中 (p) 是质数，(k) 是正整数。
显然我们可以将这样的函数拆成两个，即对于 (x=p^k)，把 (f(x)) 拆成 (f_1(x)=x) 和 (f_2(x)=x^2)，分别计算它们，求和后就是原来的 (f)。
在接下来的讨论中，我们设 (f'(x)=x^k)。分别将 (k=1,2) 带入并求和，就得到了答案。
观察到 (f') 是完全积性函数，且在质数处和 (f_k) 取值相同。
来做一些奇怪的工作
设 (mathbb P) 表示质数集， (P_i) 表示第 (i) 个质数。特别地，设 (P_0=0)。设 (L(i)) 表示 (i) 的最小质因子。考虑有这样一个数组 (g)，它满足以下定义：
其中 (j leq Q)，(Q) 是最大的满足 (P_Q leq sqrt n) 的正整数。
即，在 (1) 到 (n) 的所有 (i) 中，对所有满足 (i) 是质数或 (i) 的最小质因子大于 (P_j) 的 (i) 的 (f') 值求和。
考虑 (g(n,j)) 的转移。如果 (n<P_j^2)，那么 (n) 的最小质因子必然不会为 (P_j)，所以这个限制不会对贡献产生影响，于是有 (g(n,j)=g(n,j-1))。
如果 (n geq P_j^2) 怎么办呢？这种情况下，(g(n,j)) 会在 (g(n,j-1)) 的基础上减少。观察到 (P_{j-1},P_j) 是连续的两个质数，于是从贡献中消失的 (i) 的最小质因子一定是 (P_j)。这一部分 (i) 的贡献即是 (f'(P_j) imes( g(left lfloor dfrac{n}{P_j} ight floor,j-1)-g(P_{j-1},j-1)))。为什么呢？观察式子开头的 (f'(P_j))，这是因为 (f') 是完全积性函数，且从贡献中消失的 (i) 的最小质因子一定是 (P_j)，我们可以将这一部分的贡献提出。(g(left lfloor dfrac{n}{P_j} ight floor,j-1)) 中的 (i)（姑且记作 (i')）满足 (L(i')>P_{j-1}) 或 (i' in mathbb P)。
对于 (L(i')>P_{j-1}) 这一部分，(i') 的最小质因子大于 (P_{j-1})，那么 (i' imes P_j)（为什么要乘上 (P_j)？因为我们将它提到外面去了）的最小质因子只可能为 (P_j)，(i= i' imes P_j) 正是贡献中应该消失的 (i)。
为什么要加上 (g(P_{j-1},j-1)) 呢？这对应着第二部分 (i' in mathbb P)。对于质数 (i')，仅当 (i' imes P_j) 的最小质因子为 (P_j) 时才被删掉。观察到对于小于 (P_j) 的质数 (i')，我们不应该删去 (i =i' imes P_j) 的贡献。同时，观察到 (i') 要小于等于 (left lfloor dfrac{n}{P_j} ight floor) 才会错误地删掉贡献。因为 (P_j) 一定不大于 (sqrt n)，所以有 (P_j leq left lfloor dfrac{n}{P_j} ight floor)，这个限制可以忽略。我们要加回来的就是小于 (P_j) 的质数的贡献，即前 (j-1) 个质数的 (f') 值之和。由 (g) 的定义可知，这和 (g(P_{j-1},j-1)) 相等。
于是我们得到了 (g) 的递推式：
求解答案
设 (S(n,x)) 表示 (sum limits_{1 leq i leq n} [L(i)> P_x] f_k(i))，即 (1 sim n) 中所有最小质因子大于 (P_x) 的数的 (f_k) 值之和，注意是 (f_k) 而非 (f')。
分类讨论：
-
如果符合条件的 (i) 是质数
那么这一部分 (i) 之和为 (g(n,Q)- sum limits_{i=1}^{x} f_k(P_i))。后面的和式可以预处理的时候做一个前缀和算出来。
-
如果符合条件的 (i) 是合数
枚举最小质因子 (p_j^e)。这一部分的贡献即为
[sum limits_{p_j^e leq n and j > i }f_k(p^e) imes (S(left lfloor dfrac{n}{p_j^e} ight floor,j)+[e>1]) ](+[e>1]) 是因为 (p_j^e) 在 (e>1) 的时候本身就是合数，应该被算进去。而 (e=1) 的时候就是质数，不在这个分类里面。
综合一下，我们可以得到
这样问题就解决啦，答案就是 (S(n,0))。值得一提的是，我们没有计算 (f(1)) 的值，因此还要加上 (1)。
时间复杂度据说是 (O(dfrac{n^{frac34}}{log n}))，大概是 (1000 ext{ms}) 跑 (n=10^{10})。
# include <bits/stdc++.h>
const int N=1000010,INF=0x3f3f3f3f,MOD=1e9+7,INV6=166666668,INV2=500000004;
typedef long long ll;
ll prime[N],psum;
bool vis[N];
ll g1[N],g2[N],prsum1[N],prsum2[N],w[N],wtot,idx1[N],idx2[N];
ll MAXN;
ll n;
inline int read(void){
int res,f=1;
char c;
while((c=getchar())<'0'||c>'9')
if(c=='-')f=-1;
res=c-48;
while((c=getchar())>='0'&&c<='9')
res=res*10+c-48;
return res*f;
}
inline void init(void){ // 预处理
for(int i=2;i<=MAXN;++i){
if(!vis[i]){
prime[++psum]=i;
prsum1[psum]=(prsum1[psum-1]+i)%MOD;
prsum2[psum]=(prsum2[psum-1]+1ll*i*i)%MOD;
}
for(int j=1;i*prime[j]<=MAXN&&j<=psum;++j){
vis[i*prime[j]]=true;
if(i%prime[j]==0){
break;
}
}
}
return;
}
inline ll getsum1(ll x){ // 记得取模
x%=MOD;
return x*(x+1)%MOD*INV2%MOD;
}
inline ll getsum2(ll x){
x%=MOD;
return x*(x+1)%MOD*(2ll*x+1)%MOD*INV6%MOD;
}
inline int getid(ll x){ // trick: n/x 只有 O(sqrt n) 种取值，于是我们开大小为 2sqrt n 的数组即可，idx1[x] 存储 x <= sqrt n 时 n/x 的取值，idx2[n/x] 存储 x > sqrt n 时 n/x 的取值
return (x<=MAXN)?idx1[x]:idx2[n/x];
}
ll S(ll x,ll i){ // 递归计算 S
if(prime[i]>=x){
return 0;
}
int nowid=getid(x);
ll res=(((g2[nowid]-g1[nowid])-(prsum2[i]-prsum1[i]))%MOD+MOD)%MOD;
for(int j=i+1;j<=psum&&prime[j]*prime[j]<=x;++j){
for(ll k=1,nowp=prime[j];nowp<=x;++k,nowp=nowp*prime[j]){
ll truenowp=nowp%MOD;
res=(res+(truenowp-1)*truenowp%MOD*(S(x/nowp,j)+(k>1)))%MOD;
}
}
return res;
}
int main(void){
scanf("%lld",&n);
MAXN=1ll*sqrt(n);
init();
for(ll l=1,r;l<=n;l=r+1){
ll val=n/l;
r=n/val,w[++wtot]=val;
g1[wtot]=getsum1(val)-1,g2[wtot]=getsum2(val)-1; // 预处理出 g(...,0)
if(val<=MAXN){
idx1[val]=wtot;
}else{
idx2[n/val]=wtot;
}
}
// g1,g2 分别代表 k=1,2 时 g 数组的值
for(int i=1;i<=psum;++i){ // 滚动计算 g
for(int j=1;prime[i]*prime[i]<=w[j]&&j<=wtot;++j){
int nowid=getid(w[j]/prime[i]);
g1[j]=(g1[j]-prime[i]*(g1[nowid]-prsum1[i-1])%MOD+MOD)%MOD;
g2[j]=(g2[j]-prime[i]*prime[i]%MOD*(g2[nowid]-prsum2[i-1])%MOD+MOD)%MOD;
}
}
printf("%lld",(S(n,0)+1ll)%MOD);
return 0;
}