BZOJ3512 DZY Loves Math IV <莫比乌斯反演+杜教筛>

Problem

DZY Loves Math IV

$\mathrm{Time\;Limit:\;15\;Sec}$
$\mathrm{Memory\;Limit:\;128\;MB}$

Description

给定$n,m$，求$\sum_{i=1}^{n}\sum_{j=1}^{m}\varphi(ij)$模$10^9+7$的值。

Input

仅一行，两个整数$n,m$。

Output

仅一行答案。

Sample Input

100000 1000000000

Sample Output

857275582

HINT

$1\le n\le10^5$，$1\le m\le10^9$。

Source

By Jc
鸣谢杜教

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

Solution

杜教筛经典题。

以下解析部分转载自blutrex的题解。

由于$n$很小，考虑枚举$i$，则转化为求$S(n,m)=\sum_{i=1}^{m}\varphi(ni)$

• $n=1$时，为$\varphi$的前缀和，用杜教筛在$O(n^{\frac{2}{3}})$内求解。
• $n>1$时，若$n$无平方因子，则有

$$\varphi(ni)=\varphi(i)\sum_{d|n,i}\varphi(\frac{n}{d})$$

那么

$$\begin{aligned} S(n,m) &=\sum_{i=1}^{m}\varphi(i)\sum_{d|n,i}\varphi(\frac{n}{d})\\ &=\sum_{d|n}\varphi(\frac{n}{d})\sum_{d|i}^{m}\varphi(i)\\ &=\sum_{d|n}\varphi(\frac{n}{d})\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}\varphi(id)\\ &=\sum_{d|n}\varphi(\frac{n}{d})S(d,\lfloor\frac{m}{d}\rfloor)\\ \end{aligned}$$

如此，我们只需要一开始将$n$转化为一个无平方因子的数即可。令$k$为$n$所有质因数乘积，则有$S(n,m)=\frac{n}{k}S(k,m)$。

最终答案为$\sum_{i=1}^{n}S(i,m)$，总时间复杂度为$O(n^\frac{2}{3}+n(\sqrt{n}+\sqrt{m}))$。

Code

#include <bits/stdc++.h>
#define MAX_N 2000000
#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');
}
int cnt, pri[MAX_N+5];
int phi[MAX_N+5], mi[MAX_N+5], s[MAX_N+5];
bool NotPri[MAX_N+5]; map <lnt, lnt> ex, pr;
void init() {
    mi[1] = phi[1] = s[1] = 1;
    for (int i = 2; i <= MAX_N; i++) {
        if (!NotPri[i]) pri[cnt++] = i, phi[i] = i-1, mi[i] = i;
        for (int j = 0; j < cnt; j++) {
            if (i > MAX_N/pri[j]) break; NotPri[i*pri[j]] = true;
            if (i%pri[j]) phi[i*pri[j]] = phi[i]*(pri[j]-1), mi[i*pri[j]] = mi[i]*pri[j];
            else {phi[i*pri[j]] = phi[i]*pri[j], mi[i*pri[j]] = mi[i]; break;}
        }
        s[i] = (s[i-1]+phi[i])%P;
    }
}
lnt S(int n, int m) {
    if (m <= 1) return phi[n*m];
    if (n == 1) {
        if (m <= MAX_N) return s[m];
        lnt &ret = pr[m]; if (ret) return ret%P;
        ret = 1LL*m*(m+1)/2;
        for (int l = 2, r; l <= m; l = r+1)
            r = m/(m/l), ret -= (r-l+1)*S(1, m/l);
        return ret%P;
    }
    lnt &ret = ex[1LL*n*P+m]; if (ret) return ret;
    for (int i = 1; i*i <= n; i++) if (n%i == 0) {
        ret = (ret+1LL*phi[n/i]*S(i, m/i)%P)%P;
        if (i*i != n) ret = (ret+phi[i]*S(n/i, m/(n/i)))%P;
    }
    return ret;
}
int main() {
    int n, m; lnt ans = 0;
    read(n), read(m), init();
    for (int i = 1; i <= n; i++)
        ans = (ans+i/mi[i]*S(mi[i], m)%P)%P;
    return printf("%lld\n", ans), 0;
}