bzoj3512 DZY Loves Math IV

wennitao

于 2019-05-20 13:59:55 发布

阅读量164

点赞数

分类专栏：数论文章标签：数论莫比乌斯反演杜教筛

本文链接：https://blog.csdn.net/wennitao222/article/details/90373632

版权

数论专栏收录该内容

16 篇文章 0 订阅

订阅专栏

题面：bzoj3512

题意：求 $\sum _{i=1} ^n \sum _{j=1} ^m \varphi (ij)$ ， $\le n \le 10^5$ ， $\le m \le 10^9$

题解： $n$ 范围比较小，枚举 $n$ 。

即求 $\sum _{i=1} ^m \varphi (ni)$

设 $\prod p_i ^{a_i}$ ，令 $\prod p_i ^{a_i - 1}$ ， $\prod p_i$ ， $p q = n$

$\begin {aligned} S (n, m) &= p \sum _{i=1} ^m \varphi (qi) \\ &= p \sum _{i = 1} ^m \varphi(q) \varphi (\frac i {gcd (i, q)}) gcd (i, q) \\ &= p \sum _{i=1} ^m \varphi (i) \varphi (\frac q {gcd (i, q)}) \sum _{d \mid i, d\mid q} \varphi (d) \\ &= p \sum _{i=1} ^m \varphi (i) \sum _{d \mid i, d \mid q} \varphi (\frac q d) \\ &= p \sum _{d \mid q} \varphi (\frac q d) \sum _{i = 1} ^{\frac m d} \varphi (id) \\ &= p \sum _{d \mid q} \varphi (\frac q d) S (d, \frac m d) \end {aligned}$

#include <bits/stdc++.h>
using namespace std;
typedef long long ll ;
const int maxn = 2e6 + 10, mod = 1e9 + 7 ;
int mn[maxn], prime[maxn], tot ;
ll phi[maxn] ;
bool vis[maxn] ;
map<int, ll> p, s[100010];
void init () {
    phi[1] = 1 ;
    for (int i = 2; i < maxn; i ++) {
        if (!vis[i]) prime[++ tot] = i, phi[i] = i - 1, mn[i] = i ;
        for (int j = 1; j <= tot && i * prime[j] < maxn; j ++) {
            vis[i * prime[j]] = 1; mn[i * prime[j]] = prime[j] ;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j] ;
                break ;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1) ;
        }
    }
    for (int i = 1; i < maxn; i ++) phi[i] = (phi[i] + phi[i - 1]) % mod ;
}
ll cal (int n) {
    if (n < maxn) return phi[n] ;
    if (p[n]) return p[n] ;
    ll res = 1ll * n * (n + 1) / 2 % mod ;
    for (int i = 2, j; i <= n; i = j + 1) {
        j = n / (n / i) ;
        res = (res - 1ll * (j - i + 1) * cal (n / i) % mod) % mod ;
    }
    return p[n] = (res + mod) % mod ;
}
ll solve (int n, int m) {
    if (!m) return 0 ;
    if (m == 1) return (cal (n) - cal (n - 1) + mod) % mod ;
    if (n == 1) return cal (m) ;
    if (s[n][m]) return s[n][m] ;
    ll res = 0; int p = 1, q = 1, tmp = n ;
    vector<int> fac ;
    while (tmp > 1) {
        int x = mn[tmp]; q *= x; tmp /= x; fac.push_back (x) ;
        while (tmp % x == 0) tmp /= x, p *= x ;
    }
    int sz = fac.size() ;
    for (int i = 0; i < 1 << sz; i ++) {
        int d = 1 ;
        for (int j = 0; j < sz; j ++)
            if (i & (1 << j)) d *= fac[j] ;
        res = (res + 1ll * (cal (q / d) - cal (q / d - 1) + mod) * solve (d, m / d) % mod) % mod ;
    }
    res = res * p % mod ;
    return s[n][m] = (res + mod) % mod ;
}
int main() {
    init () ;
    int n, m ;
    scanf("%d%d", &n, &m) ;
    ll ans = 0 ;
    for (int i = 1; i <= n; i ++) ans = (ans + solve (i, m)) % mod ;
    printf("%lld\n", (ans + mod) % mod) ;
    return 0 ;
}