BZOJ 4176. Lucas的数论

 

求 $$\sum_{i=1}^n\sum_{j=1}^md(ij)$$
$$d(ij)=\sum_{x|i}\sum_{y|j}[(x,y)=1]=\sum_{x|i}\sum_{y|j}\sum_{p|(x,y)}\mu(p)$$
$$=\sum_p \mu(p)\sum_{x|i}[d|x]\sum_{y|j}[d|y]=\sum_{p|i,p|j}\mu(p)d(\frac{i}{p})d(\frac{j}{p})$$
代回原式得
$$\sum_{i=1}^n\sum_{j=1}^m\sum_{p|i,p|j}\mu(p)d(\frac{i}{p})d(\frac{j}{p})=\sum_{p}\mu(p)\sum_{p|i}d(\frac{i}{p})\sum_{p|j}d(\frac{j}{p})$$
$$=\sum_{p}\mu(p)\sum_{i=1}^{\lfloor \frac{n}{p} \rfloor}d(i)\sum_{j=1}^{\lfloor \frac{m}{p} \rfloor}d(j)$$
令 $s(n)=\sum_{i=1}^n d(i)=\sum_{i=1}^n\sum_{d=1}^i[d|i]=\sum_{d=1}^n \lfloor \dfrac{n}{d} \rfloor$
原式为$$\sum_{p}\mu(p)s(\lfloor \frac{n}{p} \rfloor)s(\lfloor \frac{m}{p} \rfloor)$$

$\mu$ 的前缀和用杜教筛，$s$ 的前缀和用整除分块

#include <bits/stdc++.h>

const int N = 5e4 + 7;
const int MOD = 1000000007;
int mu[N], prime[N], prin, d[N], t[N];
bool vis[N];

inline void M(int &ans) {
    if (ans >= MOD) ans -= MOD;
    if (ans < 0) ans += MOD;
}

void init(int n) {
    mu[1] = d[1] = t[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) {
            prime[++prin] = i;
            d[i] = 2;
            t[i] = 1;
            mu[i] = -1;
        }
        for (int j = 1; j <= prin && i * prime[j] < N; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                t[i * prime[j]] = t[i] + 1;
                mu[i * prime[j]] = 0;
                d[i * prime[j]] = d[i] / (t[i] + 1) * (t[i] + 2);
                break;
            }
            d[i * prime[j]] = d[i] * 2;
            t[i * prime[j]] = 1;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= n; i++)
        M(mu[i] += mu[i - 1]), M(d[i] += d[i - 1]);
}

std::unordered_map<int, int> muu, dd;


inline int Mu(int n) {
    if (n < N) return mu[n];
    if (muu.count(n)) return muu[n];
    int ans = 1;
    for (int i = 2, j; i <= n; i = j + 1) {
        j = n / (n / i);
        M(ans -= 1LL * (j - i + 1) * Mu(n / i) % MOD);
    }
    return muu[n] = ans;
}

inline int D(int n) {
    if (n < N) return d[n];
    if (dd.count(n)) return dd[n];
    int ans = 0;
    for (int i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        M(ans += 1LL * (j - i + 1) * (n / i) % MOD);
    }
    return dd[n] = ans;
}

int solve(int n) {
    int ans = 0;
    for (int i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        M(ans += 1LL * (Mu(j) - Mu(i - 1) + MOD) * D(n / i) % MOD * D(n / i) % MOD);
    }
    return ans;
}

int main() {
    init(N - 1);
    int n;
    scanf("%d", &n);
    printf("%d\n", solve(n));
    return 0;
}