2021 CCPC威海 I. Distance （min25筛）

题目链接
题意是给出n的哈斯图，求${\sum }_i{\sum }_{j}dis(i,j)$，$dis(i,j)$表示$i$和$j$的最短路径长度。（哈斯图为，如果$x|y$则有一条$y-x$的边，权值为$y/x$）

解题思路：
官方题解的做法非常的简约和优雅，这里记录一下推了很久常数又大的纯min25筛做法。。
我们让$f(i)$表示$i$的所有 质因数$\ast$指数 的和。
例如$f(12)=2\ast 2+3,f(48)=2\ast 4+3$
显然$f(i)$表示$i$到1的最短距离
那么有
$$ans=\sum _{i=1}^n\sum _{j=1}^{n}f(i/gcd(i,j))+f(j/gcd(i,j))$$
也就是说，$i$到$j$的最短距离 = $i$到$gcd$的最短距离+$gcd$到j的最短距离

由$f$的定义有
$$f(i/gcd)=f(i)-f(gcd)$$

所以$$ans=\sum _{i=1}^n\sum _{j=1}^n(f(i)+f(j))-2\sum _{i=1}^n\sum _{j=1}^{n}f(gcd(i,j))$$

假设我们可以快速计算$s(n)={\sum }_{i=1}^{n}f(i)$

那么
$$ans=2n\ast s(n)-2\sum _{d=1}^n\sum _{i=1}^{n/d}\sum _{j=1}^{n/d}f(d)\ast [gcd(i,j)==1]$$
又因为
$$\sum _{d=1}^n\sum _{i=1}^{n/d}\sum _{j=1}^{n/d}f(d)\ast [gcd(i,j)==1]=\sum _{d=1}^{n}f(d)\ast (2\sum _{i=1}^{n/d}\varphi (i)-1)$$
于是，如果我们可以快速算出s(n)，我们就可以筛出欧拉函数前缀和之后进行整除分块

现在问题只剩下如何快速求出s(n)

考虑$min25$筛
用$F(k,n)$表示 ${\sum }_{lnp(i)\ge p_k||i\in Prime}f(i)$
用$h(k,n)$表示 ${\sum }_{lnp(i)\ge p_k||i\in Prime}1$，也就是满足条件的n以内的数字个数
其中$lnp(i)$表示i的最小质因子，$p_k$表示第k大的质因子
那么考虑枚举$p_k$的次数进行转移
$$F(k,n)=F(k+1,n)+\sum _{c\ge 1\&\&p_k^{c+1}\le n}c\ast p_k(h(k+1,n/p_k^c)-h(k+1,p_k))+F(k+1,n/p_k^c)-F(k+1,p^k)+(c+1)p_k$$
$$h(k,n)=h(k+1,n)+\sum _{c\ge 1\&\&p_k^{c+1}\le n}(h(k+1,n/p_k^c)-h(k+1,p_k^c))+1$$
记得设置初始值$F(K,n)，h(K,n)$，表示只有质数的时候对应的值。
顺便也把欧拉函数前缀和筛一下，然后进行整除分块。
但是常数很大，时间花费主要是在min25的$sol$这部分。在cf上得交c++17(64)才能过
这篇blog单纯记录一下这种转移方法方便以后回想QAQ

#include<bits/stdc++.h>
#define fors(i,a,b) for(int i = a; i < b; ++i)
#define ll long long
using namespace std;
const int mod = 1e9+7;
const int inv2 = (mod+1)/2;
namespace min25{
    const int maxn = 326232;
    const int N = maxn;
    int p[maxn], cnt = 0, H[maxn];
    void INIT(){
        fors(i,2,maxn){
            if(!H[i]) p[++cnt] = i;
            for(int j = 1; j <= cnt && p[j]*i < maxn; ++j){
                H[i*p[j]] = 1; 
                if(i%p[j] == 0)  break;
            }
        }
    }
    ll n, B;
    int g1[N*2], g2[N*2];
    ll w[N*2];
    int id1[N], id2[N], tot;
    int id(ll x){return x <= B?id1[x]:id2[n/x];}
    void init(ll _n){
        n = _n; B = sqrt(n); tot = 0;
        for(ll i = 1,j; i <= n; i = j+1){
            j = n/(n/i); w[++tot] = n/i; 
            if(w[tot] <= B) id1[w[tot]] = tot;
            else id2[n/w[tot]] = tot;
            g1[tot] = (w[tot]-1)%mod;
            g2[tot] = (w[tot]-1)%mod*((w[tot]+2)%mod)%mod * inv2 %mod;
        }
        for(int i = 1; i <= cnt; ++i){
            for(int j = 1; j <= tot && (ll)p[i]*p[i] <= w[j]; ++j){
                int k1 = id(w[j]/p[i]), k2 = id(p[i-1]);
                g1[j] = (g1[j] - ((ll)g1[k1]-g1[k2]))%mod;
                g2[j] = (g2[j] - (ll)p[i]*(g2[k1] - g2[k2]))%mod;
            }
        }
        return;
    }
    int f[N*2], h[N*2], sp[N*2];
    void sol(){
        for(int i = 1; i <= tot; ++i) f[i] = g2[i], h[i] = g1[i], sp[i] = (g2[i]-g1[i])%mod;
        for(int i = cnt; i >= 1; --i){
            for(int j = 1; j <= tot && (ll)p[i]*p[i] <= w[j]; ++j){
                ll cur = p[i];
                ll phi = p[i]-1;
                for(int e = 1; cur*p[i] <= w[j]; ++e, cur *= p[i], phi *= p[i]){
                    f[j] = (f[j] + f[ id(w[j]/cur) ] - (ll)f[id(p[i])] + (ll)(e+1)*p[i] + (ll)e*p[i]%mod* (h[id(w[j]/cur)] - h[id(p[i])]) )%mod;
                    sp[j] = ( sp[j] + (ll)phi*( sp[id(w[j]/cur)] - sp[id(p[i])] )%mod + (ll)phi*p[i] )%mod;
                    h[j] = (h[j] + (ll)h[ id(w[j]/cur)]-h[id(p[i])]+1)%mod;
                }
            }
        }
    }
    inline ll get_f(ll n){
        if(n < 2) return 0;
        return f[id(n)];
    }
    inline ll get_phi(ll n){
        if(n < 2) return n == 1;
        return sp[id(n)]+1;
    }
}
int main(){
    ll n; cin>>n;
    min25::INIT();
    min25::init(n);
    min25::sol();
    ll ans = (2LL*n+2)%mod*min25::get_f(n)%mod;
    ll l = 1;
    while(l <= n){
        ll k = n/l;
        ll r = n/k;
        ll sf = min25::get_f(r)-min25::get_f(l-1);
        sf %= mod;
        ans = (ans-2LL*sf%mod*(2LL*min25::get_phi(k)))%mod;
        l = r+1;
    }
    ans = (ans + mod)%mod;
    cout<<ans<<endl;
    return 0;
}