洛谷 P1829 [国家集训队]Crash的数字表格（莫比乌斯反演）

传送门

---

题目大意

求${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$

常规解法

下面表述默认$n\le m$

$f(n,m)={\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)={\sum }_{i=1}^n{\sum }_{j=1}^m\frac{ij}{gcd(i,j)}$

有一个常用的技巧，首先枚举$gcd(i,j)$：$f(n,m)={\sum }_{d=1}^n{\sum }_i^n{\sum }_j^m\frac{ij}{d}[gcd(i,j)=d]$

按照套路，$[gcd(i,j)=d]\iff [gcd(\lfloor \frac{i}{d}\rfloor ,\lfloor \frac{j}{d}\rfloor )=1]$，从而转化为内层的两个求和同除以$d$，又因为外面乘上了$ij$，那么需要再多乘上$d^2$，得：$f(n,m)={\sum }_{d=1}^n{\sum }_i^{\lfloor \frac{n}{d}\rfloor }{\sum }_j^{\lfloor \frac{m}{d}\rfloor }i\ast j\ast d\ast [gcd(i,j)=1]$

然后对于$[gcd(i,j)=1]$，依然可以利用莫比乌斯变换替换：$f(n,m)={\sum }_{d=1}^n{\sum }_i^{\lfloor \frac{n}{d}\rfloor }{\sum }_j^{\lfloor \frac{m}{d}\rfloor }i\ast j\ast d{\sum }_{k|(i,j)}\mu (k)$

枚举$k$，类似上面那样令$i=x\ast k,j=y\ast k$：$f(n,m)={\sum }_{d=1}^{n}d{\sum }_{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu (k)\ast k^2{\sum }_{x=1}^{\lfloor \frac{n}{dk}\rfloor }x{\sum }_{y=1}^{\lfloor \frac{m}{dk}\rfloor }y$

两次数论分块时间复杂度$O(n)$，${\sum }_d^{n}d,{\sum }_i^{n}i,{\sum }_j^{m}j$都可以使用等差数列通项公式

这里的两次数论分块看了很久，一定要细细体会

代码

//
// Created by Happig on 2020/9/24
//
#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 = 20101009;
const int maxn = 1e7 + 10;

vector<int> prime;
int mu[maxn];
ll sum[maxn];
bitset<maxn> vis;

void init() {
    vis.reset(), prime.clear();
    mu[1] = sum[1] = 1;
    for (int i = 2; i < maxn; i++) {
        if (!vis[i]) {
            prime.push_back(i);
            mu[i] = -1;
        }
        for (int j = 0; j < prime.size() && i * prime[j] < maxn; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j]) {
                mu[i * prime[j]] = -mu[i];
            } else {
                mu[i * prime[j]] = 0;
                break;
            }
        }
        sum[i] = (sum[i - 1] + (1LL * mu[i] * i * i) % Mod) % Mod;
    }
}

ll g(ll x, ll y) {
    return (x * (1 + x) / 2) % Mod * ((y * (1 + y) / 2) % Mod) % Mod;
}

ll cal(ll n, ll m) {
    ll up = min(n, m), ans = 0;
    for (ll l = 1, r; l <= up; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        if (r > up) r = up;
        ans += (sum[r] - sum[l - 1] + Mod) % Mod * g(n / l, m / l) % Mod;
    }
    return ans;
}

ll solve(ll n, ll m) {
    ll up = min(n, m), ans = 0;
    for (ll l = 1, r; l <= up; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        if (r > up) r = up;
        ans += (r - l + 1) * (l + r) / 2 % Mod * cal(n / l, m / l) % Mod;
        ans %= 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);
    int n, m;
    init();
    cin >> n >> m;
    cout << solve(n, m) << ENDL;
    return 0;
}

继续优化

上式还可以继续化简优化：

设$T=dk,f(x)={\sum }_{i=1}^{x}i$，那么可以得到：$f(n,m)={\sum }_{d=1}^{n}d{\sum }_{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu (k)\ast k^{2}f(\lfloor \frac{n}{T}\rfloor )f(\lfloor \frac{m}{T}\rfloor )$

然后仍使用上面枚举$gcd(i,j)$的技巧，枚举$T$（这里比较难懂，因为$T$是$d$的倍数，然后$T$可以取遍$[1,n]$，故这样转化）：
$$\begin{gathered} f(n,m)={\sum }_{T=1}^{n}f(\lfloor \frac{n}{T}\rfloor )f(\lfloor \frac{m}{T}\rfloor ){\sum }_{d|T}d\ast \mu (\frac{T}{d})\ast (\frac{T}{d}{)}^2 \\ ={\sum }_{T=1}^{n}f(\lfloor \frac{n}{T}\rfloor )f(\lfloor \frac{m}{T}\rfloor ){\sum }_{d|T}d\ast \mu (d)\ast T \end{gathered}$$

上面的${\sum }_{d|T}d\ast \mu (\frac{T}{d})\ast (\frac{T}{d}{)}^2\iff {\sum }_{d|T}d\ast \mu (d)\ast T$比较巧妙，原理是：

• 当$d$取遍$T$的所有因数时，因为$d|T$，那么$\mu (\frac{T}{d})=\mu (d)$
• 当$d$取遍$T$的所有因数时，那么$d\ast (\frac{T}{d}{)}^2=\frac{T^2}{d}=T\ast \frac{T}{d}=T\ast d$

设$F(T)={\sum }_{d|T}d\ast \mu (d)$，考虑$a\perp b$：$F(a)={\sum }_{d|a}d\ast \mu (d),F(b)={\sum }_{d|b}d\ast \mu (d)$

在我们枚举$d$时，它既可以是$a$的约数也可以是$b$的约数，根据唯一分解定理不难发现$F(ab)=F(a)\ast F(b)$，故$F(n)$是一个积性函数，可以线性筛。可以知道$F(1)=1,F(p)=1-p,p\in primes$；考虑$a\%p==0$的情况，那么$a\ast p$中分解出的因数含$p^{\alpha +\beta }$项的值一定为$0$，因为其幂次一定大于$1$，而大于$1$之后莫比乌斯函数的值即为$0$，因此可以得出此时$F(a\ast p)=F(a),a\%==0$，现在就可以欧拉筛预处理出$F(n)$了

此时再观察公式${\sum }_{T=1}^{n}f(\lfloor \frac{n}{T}\rfloor )f(\lfloor \frac{m}{T}\rfloor )F(T)\ast T$，那么$F(T)\ast T$也可以作为前缀和预处理出来了，那么算法时间复杂度优化了很多

代码

//
// Created by Happig on 2020/9/24
//
#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 = 20101009;
const int maxn = 1e7 + 10;

ll sum[maxn], F[maxn];
bitset<maxn> vis;
vector<int> prime;

void init() {
    sum[1] = F[1] = 1;
    vis.reset(), prime.clear();
    for (int i = 2; i < maxn; i++) {
        if (!vis[i]) {
            prime.push_back(i);
            F[i] = (1 - i + Mod);
        }
        for (int j = 0; j < prime.size() && i * prime[j] < maxn; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j]) {
                F[i * prime[j]] = (F[i] * F[prime[j]]) % Mod;
            } else {
                F[i * prime[j]] = F[i];
                break;
            }
        }
        sum[i] = (sum[i - 1] + F[i] * i % Mod) % Mod;
    }
}

ll f(ll x) {
    return x * (x + 1) / 2 % Mod;
}

ll solve(ll n, ll m) {
    ll up = min(n, m), ans = 0;
    for (ll l = 1, r; l <= up; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        if (r > up) r = up;
        ans += (sum[r] - sum[l - 1] + Mod) % Mod * f(n / l) % Mod * f(m / l) % Mod;
        ans %= 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);
    int n, m;
    init();
    cin >> n >> m;
    cout << solve(n, m) << ENDL;
    return 0;
}