LibreOJ - 6682 梦中的数论（min_25筛）

传送门

---

题目大意

求${\sum }_{i=1}^n{\sum }_{j=1}^n{\sum }_{k=1}^n[(j|i)\&\&(j+k)|i]$，即求$j$整除$i$且$j+k$整除$i$的个数的前缀和

解题思路

看到$j|i$不难想到$j$是$i$的所有因数，而$(j+k)|i$实际上上的个数即为大于$j$的$i$的因子的个数。那么实际上就是对所有因子选出不重复的对数的个数，设$d(i)$为$i$的因子个数，那么原式可化为：

${\sum }_{i=1}^{n}C(d(i),2)={\sum }_{i=1}^n\frac{d(i{)}^2-d(i)}{2}$

对于${\sum }_{i=1}^{n}d(i)$，实际上就是${\sum }_{i=1}^n\lfloor \frac{n}{i}\rfloor$，整除分块即可；对于${\sum }_{i=1}^{n}d(i{)}^2$，因为$d(i)$为积性函数，那么$f(i)=d(i{)}^2$也为积性函数，又因为$f(p)=4,f(p^k)=(k+1{)}^2$，均可以快速计算，那么考虑亚线性筛法$min\_25$筛，直接按模板套一下就好了

超详细的min_25筛

//
// Created by Happig on 2020/10/6
//
#include <bits/stdc++.h>
#include <unordered_map>
#include <unordered_set>

using namespace std;
#define fi first
#define se second
#define pb push_back
#define ins insert
#define Vector Point
#define ENDL "\n"
#define lowbit(x) (x&(-x))
#define mkp(x, y) make_pair(x,y)
#define mem(a, x) memset(a,x,sizeof a);
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef pair<double, double> pdd;
const double eps = 1e-8;
const double pi = acos(-1.0);
const int inf = 0x3f3f3f3f;
const double dinf = 1e300;
const ll INF = 1e18;
const int Mod = 998244353;
const int maxn = 2e5 + 10;

ll prime[maxn];
ll w[maxn], id1[maxn], id2[maxn], g[maxn];
bitset<maxn> vis;
int m, cnt, tot;
ll n;

ll qkp(ll x, ll n, ll p) {
    ll ans = 1;
    x %= p;
    while (n) {
        if (n & 1) ans = ans * x % p;
        x = x * x % p;
        n >>= 1;
    }
    return ans;
}

ll getInv(ll a, ll p) {
    return qkp(a, p - 2, p);
}

ll cal() {
    ll ans = 0;
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        if (r > n) r = n;
        ans += (r - l + 1) * (n / l) % Mod;
        ans %= Mod;
    }
    return ans;
}

ll getID(ll x) {
    return x <= m ? id1[x] : id2[n / x];
}

void init() {
    vis.reset();
    vis[1] = 1;
    m = sqrt(n + 0.5), cnt = tot = 0;
    for (int i = 2; i <= m; i++) {
        if (!vis[i]) prime[++cnt] = i;
        for (int j = 1; j <= cnt && 1LL * i * prime[j] <= m; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }

    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        w[++tot] = n / l;
        g[tot] = w[tot] - 1;
        if (w[tot] <= m) id1[w[tot]] = tot;
        else id2[r] = tot;
    }

    for (int i = 1; i <= cnt; i++) {
        for (int j = 1; j <= tot && prime[i] * prime[i] <= w[j]; j++) {
            ll k = getID(w[j] / prime[i]);
            g[j] -= g[k] - i + 1;
        }
    }
}

ll S(ll x, int y) {
    if (x <= 1 || prime[y] > x) return 0;
    ll k = getID(x);
    ll ans = (g[k] - y + 1) * 4 % Mod;
    for (int i = y; i <= cnt && prime[i] * prime[i] <= x; i++) {
        ll pe = prime[i];
        for (int e = 1; pe * prime[i] <= x; e++, pe *= prime[i]) {
            ans = (ans + 1LL * (e + 1) * (e + 1) % Mod * S(x / pe, i + 1) % Mod + (e + 2) * (e + 2) % Mod) % Mod;
        }
    }
    return ans;
}

int main() {
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    ll inv = getInv(2, Mod);
    cin >> n;
    init();
    ll ans1 = S(n, 1) + 1, ans2 = cal();
    //cout << ans1 << " " << ans2 << ENDL;
    cout << (ans1 - ans2 + Mod) % Mod * inv % Mod << ENDL;
    return 0;
}