BZOJ 4407 - 于神之怒加强版

#数学 #莫比乌斯反演 #线性筛
Published on 2016-08-12

题目地址

描述

给定 k(1≤k≤5000000)k(1\le k\le 5000000)，另有 T(T≤2000)T(T\le 2000) 组 n,m(1≤n,m≤5000000)n, m(1\le n, m\le 5000000)，求：

∑i=1n∑j=1mgcd(i,j)kmod(109+7)\sum_{i = 1}^n\sum_{j = 1}^m\gcd(i, j)^k \bmod ({10}^9 + 7)

分析

直接上老一套，枚举最大公约数 dd 然后转换为可分块的形式：

如果设 f(n)=∑d∣ndkμ(Td)f(n) = \sum_{d \mid n}d^k\mu(\frac T d)，则我们只需要求 ff 的一个前缀和就可以 O(n)O(\sqrt n) 回答询问了。
发现 dkd^k 是积性函数，根据莫比乌斯反演的性质，ff 也是一个积性函数，我们可以线筛 O(n)O(n) 求出答案。
这里化一下 ff，免得线筛卡住了：

总复杂度：O(n+Tn)O(n + T\sqrt n)。

代码

//  Created by Sengxian on 8/12/16.
//  Copyright (c) 2016年 Sengxian. All rights reserved.
//  BZOJ 4407 莫比乌斯反演
# pragma GCC optimize("O3")
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
inline int ReadInt() {
    static int n, ch;
    n = 0, ch = getchar();
    while (!isdigit(ch)) ch = getchar();
    while (isdigit(ch)) n = (n << 3) + (n << 1) + ch - '0', ch = getchar();
    return n;
}

const int maxn = 5000000 + 3, modu = 1000000007;
int mu[maxn], sum[maxn], primes[maxn], num;
bool is_prime[maxn];
int K;

inline int pow_mod(ll a, int b) {
    ll res = 1;
    while (b) {
        if (b & 1) (res *= a) %= modu;
        (a *= a) %= modu;
        b >>= 1;
    }
    return res;
}

void sieve() {
    fill(is_prime, is_prime + maxn, 1);
    sum[1] = 1, num = 0;
    for (int i = 2; i < maxn; ++i) {
        if (is_prime[i]) primes[num++] = i, sum[i] = (pow_mod(i, K) - 1) % modu;
        for (int j = 0; j < num && i * primes[j] < maxn; ++j) {
            static int d;
            d = i * primes[j];
            is_prime[d] = false;
            if (i % primes[j] == 0) {
                sum[d] = (ll)sum[i] * (sum[primes[j]] + 1) % modu;
                break;
            } else sum[d] = (ll)sum[i] * sum[primes[j]] % modu;
        }
    }
    sum[0] = 0;
    for (int i = 1; i < maxn; ++i) (sum[i] += sum[i - 1]) %= modu;
}

inline int F(int n, int m) {
    if (n > m) swap(n, m);
    int ans = 0;
    for (int i = 1, last = 1; i <= n; i = last + 1) {
        last = min(n / (n / i), m / (m / i));
        (ans += (ll)(n / i) * (m / i) % modu * (sum[last] - sum[i - 1]) % modu) %= modu;
    }
    return (ans % modu + modu) % modu;
}

int main() {
    int caseNum = ReadInt();
    K = ReadInt();
    sieve();
    while (caseNum--) {
        int n = ReadInt(), m = ReadInt();
        printf("%d\n", F(n, m));
    }
    return 0;
}
BZOJ 3994 - [SDOI2015]约数个数和