BZOJ4176 Lucas的数论 <莫比乌斯反演+杜教筛>

Problem

Lucas的数论

$\mathrm{Time\;Limit:\;30\;Sec}$
$\mathrm{Memory\;Limit:\;256\;MB}$

Description

去年的Lucas非常喜欢数论题，但是一年以后的Lucas却不那么喜欢了。
在整理以前的试题时，发现了这样一道题目“求$\sum{f(i)}$,其中$1\le i\le N$”，其中$f(i)$表示$i$的约数个数。他现在长大了，题目也变难了。
求如下表达式的值：

$$\sum_{i=1}^{n}\sum_{j=1}^{n}f(ij)$$

其中$f(ij)$表示$ij$的约数个数。
他发现答案有点大，只需要输出模$1000000007$的值。

Input

第一行一个整数$n$。

Output

一行一个整数$ans$，表示答案模$1000000007$的值。

Sample Input

2

Sample Output

8

HINT

$n\le10^9$

标签：莫比乌斯反演 杜教筛

Solution

题意和BZOJ3994相同，只是将数据范围扩大了。
仍然要用相同的结论：$d(ij)=\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]$。
然后就可以推反演了：

$$\begin{aligned} \sum_{i=1}^{n}\sum_{j=1}^{n}d(ij) &=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{p|i}\sum_{q|j}[\gcd(p,q)=1]\\ &=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{p|i}\sum_{q|j}\sum_{d|p,q}\mu(d)\\ &=\sum_{p=1}^{n}\sum_{q=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{q}\rfloor}\sum_{d|p,q}\mu(d)\\ &=\sum_{d=1}^{n}\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{q=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{n}{pd}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{qd}\rfloor}\mu(d)\\ &=\sum_{d=1}^{n}\mu(d)\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{n}{pd}\rfloor}\sum_{q=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{qd}\rfloor}1\\ &=\sum_{d=1}^{n}\mu(d)(\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{n}{pd}\rfloor}1)^2\\ &=\sum_{d=1}^{n}\mu(d)(\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{pd}\rfloor)^2\\ \end{aligned}$$

对于后面的和式，可以用数论分块在$\mathcal{O}(\sqrt{n})$的时间内算出。对前面也进行数论分块，发现$d$的范围比较大，不能直接预处理$\mu$的前缀和，于是用杜教筛计算$\mu$的前缀和即可。杜教筛可参考BZOJ3944。
总复杂度…$\mathcal{O}(n^\frac{3}{4})$？

Code

#include <bits/stdc++.h>
#define MX 2500000
#define P 1000000007
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 mu[MX+5]; map <lnt, lnt> ex;
bool NotPri[MX+5]; int pri[MX+5], cnt;
void init() {
    mu[1] = 1;
    for (int i = 2; i <= MX; i++) {
        if (!NotPri[i]) pri[cnt++] = i, mu[i] = -1;
        for (int j = 0; j < cnt; j++) {
            if (i*pri[j] > MX) break;
            NotPri[i*pri[j]] = true;
            if (i%pri[j]) mu[i*pri[j]] = -mu[i];
            else {mu[i*pri[j]] = 0; break;}
        }
        mu[i] = (mu[i]+mu[i-1])%P;
    }
}
lnt sum(lnt n) {
    if (n <= MX) return mu[n];
    if (ex[n]) return ex[n]; lnt ret = 1;
    for (lnt l = 2, r; l <= n; l = r+1)
        r = n/(n/l), ret = (ret-(r-l+1)*sum(n/l)%P)%P;
    return ex[n] = ret;
}
lnt calc(lnt n) {
    lnt ret = 0;
    for (lnt l = 1, r; l <= n; l = r+1)
        r = n/(n/l), ret = (ret+(r-l+1)*(n/l)%P)%P;
    return ret*ret%P;
}
int main() {
    lnt n, ans = 0; read(n), init();
    for (lnt l = 1, r; l <= n; l = r+1)
        r = n/(n/l), ans = (ans+(sum(r)-sum(l-1))%P*calc(n/l)%P)%P;
    return printf("%lld\n", (ans+P)%P), 0;
}