BZOJ2154 Crash的数字表格 <莫比乌斯反演>

Problem

Crash的数字表格

$\mathrm{Time\;Limit:\;20\;Sec}$
$\mathrm{Memory\;Limit:\;259\;MB}$

Description

今天的数学课上，$\mathrm{Crash}$小朋友学习了最小公倍数$\mathrm{(Least\;Common\;Multiple)}$。对于两个正整数$a$和$b$，$\mathrm{LCM(a, b)}$表示能同时被$a$和$b$整除的最小正整数。例如，$\mathrm{LCM(6, 8)=24}$。回到家后，$\mathrm{Crash}$还在想着课上学的东西，为了研究最小公倍数，他画了一张$N\times M$的表格。每个格子里写了一个数字，其中第$i$行第$j$列的那个格子里写着数为$\mathrm{LCM(i, j)}$。看着这个表格，$\mathrm{Crash}$想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当$N$和$M$很大时，$\mathrm{Crash}$就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，$\mathrm{Crash}$只想知道表格里所有数的和$\mod 20101009\;\;$的值。

Input

输入的第一行包含两个正整数，分别表示$N$和$M$。

Output

输出一个正整数，表示表格中所有数的和$\mod 20101009\;\;$的值。

Sample Input

4 5

Sample Output

122

HINT

$100\%$的数据满足$N,M\le10^7$。

标签：莫比乌斯反演

Solution

$$\begin{aligned} Ans&=\sum_{x=1}^{n}\sum_{y=1}^{m}\frac{x\times y}{\gcd(x,y)}\\ &=\sum_{d=1}^{\min(n,m)}\sum_{x=1}^{n}\sum_{y=1}^{m}[\gcd(x,y)=d]\frac{x\times y}{d}\\ \end{aligned}$$

$$\begin{aligned} &Let\;f(n,m,d)=\sum_{x=1}^{n}\sum_{y=1}^{m}[\gcd(x,y)=d]\cdot x\cdot y\\ &then\;f(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor,1)=\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(x,y)=1]\cdot x\cdot y\\ &\therefore Ans=\sum_{d=1}^{\min(n,m)}\frac{f(n,m,d)}{d}=\sum_{d=1}^{\min(n,m)}f(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor,1)\times d\\ &\;\;\;\;If\;we\;can\;calculate\;f(p,q,1)\;in\;a\;rapid\;way,\\&\;\;\;\;we\;can\;calculate\;the\;answer\;rapidly.\\ \end{aligned}$$

$$\begin{aligned} f(p,q,1)&=\sum_{i=1}^{p}\sum_{j=1}^{q}\sum_{d|\gcd(i,j)}\mu(d)\cdot i\cdot j\\ &=\sum_{d=1}^{\min(p,q)}\mu(d)\times d^2\;\sum_{i=1}^{\lfloor\frac{p}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{q}{d}\rfloor}i\times j\\ &=\sum_{d=1}^{\min(p,q)}\mu(d)\times d^2\times\frac{\lfloor\frac{p}{d}\rfloor\times(\lfloor\frac{p}{d}\rfloor+1)}{2}\times\frac{\lfloor\frac{q}{d}\rfloor\times(\lfloor\frac{q}{d}\rfloor+1)}{2}\\ \end{aligned}$$

由以上推导，可见$f(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor,1)$和$\frac{\lfloor\frac{p}{d}\rfloor\times(\lfloor\frac{p}{d}\rfloor+1)}{2}\times\frac{\lfloor\frac{q}{d}\rfloor\times(\lfloor\frac{q}{d}\rfloor+1)}{2}$是可以根号分块的，在外层对$f(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor,1)$进行分块，在每个值相同的块中，对$\frac{\lfloor\frac{p}{d}\rfloor\times(\lfloor\frac{p}{d}\rfloor+1)}{2}\times\frac{\lfloor\frac{q}{d}\rfloor\times(\lfloor\frac{q}{d}\rfloor+1)}{2}$进行分块以求出$f(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor,1)$带回外层算贡献。

综上，外层复杂度为$O(\sqrt{n})$，内层复杂度为$O(\sqrt{n})$，总时间复杂度为$O(n)$。

其实可以做得更快，详见加强版BZOJ2693。

Code

#include <bits/stdc++.h>
#define MAX_N 10000000
#define MOD 20101009
using namespace std;
typedef long long lnt;
template <class T>
inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
lnt n, m, ans, s[MAX_N+5];
int cnt, pri[MAX_N+5], mu[MAX_N+5];
bool NotPri[MAX_N+5];
void getMu() {
    NotPri[1] = true, mu[1] = 1;
    for (int i = 2; i <= min(n,m); i++) {
        if (!NotPri[i]) pri[cnt++] = i, mu[i] = -1;
        for (int j = 0; j < cnt; j++) {
            if (i*pri[j] > min(n,m)) break;
            NotPri[i*pri[j]] = true;
            if (i%pri[j]) mu[i*pri[j]] = -mu[i];
            else {mu[i*pri[j]] = 0; break;}
        }
    }
    for (int i = 1; i <= min(n,m); i++) s[i] = (s[i-1]+1LL*mu[i]*i*i%MOD)%MOD;
}
lnt f(lnt p, lnt q) {
    lnt ret = 0;
    for (lnt l = 1, r; l <= min(p, q); l = r+1)
        r = min(p/(p/l), q/(q/l)), 
        (ret += (p/l*(p/l+1)/2%MOD)*(q/l*(q/l+1)/2%MOD)%MOD*(s[r]-s[l-1])%MOD) %= MOD;
    return ret;
}
int main() {
    read(n), read(m), getMu();
    for (lnt l = 1, r; l <= min(n, m); l = r+1)
        r = min(n/(n/l), m/(m/l)), (ans += (l+r)*(r-l+1)/2%MOD*f(n/l, m/l)%MOD) %= MOD;
    return printf("%lld", (ans+MOD)%MOD), 0;
}