正题
题目链接:https://www.luogu.com.cn/problem/P5325
题目大意
定义一个积性函数满足(f(p^k)=p^k(p^k-1))
求(sum_{i=1}^nf(i))
解题思路
首先我们可以把(f(p^k))是质数的情况拆成一个(2)阶的多项式(f(x)=x^2-x)。
然后就是( ext{Min25})筛了。我们先需要求出质数有关的答案，考虑把这个多项式的多个阶分开来求。
有关质数的求法，我们定义(dp)数组
（上面(Pri)表示质数集，(minp(x))表示(x)的最小质因子，(p_j)表示(leq sqrt{n})的第(j)个质因子，(k)表示多项式的某个阶，后面为了方便(g_k)写作(g)，后文中的(g(n,j))中的(n)都与输入的(n)无关，只是一个变量）
具体的说就是如果(i)是质数或者最小质因子超过(p_j)，那当某个(p_j)最接近但小于(sqrt n)时，(g(n,j))就是我们需要的答案了。
因为一个合数(x)满足(minp(x)leqsqrt x)，所以(p_j)的数量级很小，可以利用这个性质。
递推这个数组时对于每次(j)往上加相当于多加上了一些限制条件，也就是(g(n,j))相对于(g(n,j-1))我们需要减去一些多余的答案。
这些多余的答案就是最小质因数是(p_j)的数，这些数就是(p_j^k*{large(} g(frac{n}{p_j},j-1)-g(p_{j-1},j-1) {large)})。
(g(frac{n}{p_j},j-1)-g(p_{j-1},j-1) {large)})表示枚举一个乘上(p_{j})之后不会超过(n)且最小质因子超过(p_{j})的数，这些数乘上(p_j)后就是不重不漏的表示最小质因子是(p_j)的数，因为(y=x^k)是一个完全积性函数，所以可以直接乘上一个(p_j^k)。
现在我们就有递推式
快速处理(g(n,0))然后递推，可以注意到(g(p_{j-1},j-1))就是前(j-1)个质数的(k)次方和，因为(p_j)数量级只有(sqrt n)所以可以直接预处理，因为后面还要用就定义(sp_i=sum_{i=1}^nf(p_i))
还有发现无论如何每个(g(x,j))的取值都只与(lfloorfrac{n}{x} floor)有关，所以我们对于(g)数组的每一个(j)都只有(sqrt n)个连续的范围内的同一取值。所以我们可以直接整除分块来大大缩减数量级的空间和时间。
给出一个(x)如何快速得到(g(x,j))的储存地址？一个比较巧妙的方法是如果(xleq sqrt n)那么(ind1_x)表示(x)的储存位置，否则用(ind2_{frac{n}{x}})表示(x)的储存地址。
然后(j)那个维度我们只需要用到(j=cnt)（(cnt)表示(sqrt n)以内的质数数量）值，所以我们直接滚动来递推就好了。
现在我们知道了所有的(g(i,cnt))，这道题的话具体的值就是(g_2(i,cnt)-g_1(i,cnt))就可以表示前缀质数(f)的函数和了。下面简写(g(i)=g_2(i,cnt)-g_1(i,cnt))。
然后就是要求答案了，同样使用(dp)的技巧，设
然后对于(S(n,j))的答案我们可以分为质数和非质数两部分，质数部分显然就是(g(n)-sp_n)（不记得(sp)的去看前面定义）。
然后合数部分我们考虑递归下去处理就是(sum_{k>j}^{p_{k}^eleq n}f(p_k^e){large(}S(frac{n}{p_k^e},k)+[e
eq1]{large)})。(p_k^e)是枚举质因子就不多说了，那个([e
eq 1])是因为(p_k)已经作为质数被计算过了，但是后面的(p_k^e)还没有。
现在就有(S(n,j))的递推式了
递归下去做就好了，说是不知道因为啥定理所以是不用记忆化的。
一个细节就是(f(1))最好不要统计进去，最好加上就好了
据说时间复杂度是(O(frac{n^{frac{3}{4}}}{log n}))的，反正这题我没开(O2)的话(1e10)只跑了七百多毫秒。
code
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
const ll N=1e6+10,P=1e9+7;
ll n,T,cnt,tot,pri[N],sp1[N],sp2[N],ind1[N],ind2[N],g1[N],g2[N],w[N];
bool v[N];
void init(ll n){
for(ll i=2;i<=n;i++){
if(!v[i]){
pri[++cnt]=i;
sp1[cnt]=(sp1[cnt-1]+i)%P;
sp2[cnt]=(sp2[cnt-1]+i*i)%P;
}
for(ll j=1;j<=cnt&&i*pri[j]<=n;j++){
v[i*pri[j]]=1;
if(i%pri[j]==0)break;
}
}
return;
}
ll S(ll x,ll y){
if(pri[y]>=x)return 0;
ll pos=(x>T)?ind2[n/x]:ind1[x];
ll ans=g2[pos]-g1[pos]-(sp2[y]-sp1[y]);
ans=(ans%P+P)%P;
for(ll i=y+1;i<=cnt&&pri[i]*pri[i]<=x;i++){
for(ll j=1,p=pri[i];p<=x;j++,p=p*pri[i]){
ll w=p%P;
(ans+=w*(w-1)%P*(S(x/p,i)+(j!=1))%P)%=P;
}
}
return ans;
}
signed main()
{
scanf("%lld",&n);
T=sqrt(n);init(T);
ll inv2=(P+1)/2,inv6=(P+1)/6;
for(ll l=1,r;l<=n;l=r+1){
ll x=n/l;r=n/(n/l);
w[++tot]=n/l;x%=P;
g1[tot]=x*(x+1)%P*inv2%P-1;
g2[tot]=x*(x+1)%P*(2*x+1)%P*inv6%P-1;
if(n/l<=T)ind1[n/l]=tot;
else ind2[n/(n/l)]=tot;
}
for(ll i=1;i<=cnt;i++)
for(ll j=1;j<=tot&&pri[i]*pri[i]<=w[j];j++){
ll k=(w[j]/pri[i]);k=(k>T)?ind2[n/k]:ind1[k];
(g1[j]+=P-pri[i]*(g1[k]-sp1[i-1])%P)%=P;
(g2[j]+=P-pri[i]*pri[i]%P*(g2[k]-sp2[i-1])%P)%=P;
}
printf("%lld
",(S(n,0)+1)%P);
return 0;
}