P4449 于神之怒加强版

P4449 于神之怒加强版

$$
\begin {aligned}
ANS & = \sum _ {i = 1} ^ n \sum _ {j = 1} ^ m (i, j) ^ k \\
& = \sum _ {d = 1} ^ n d ^ k \sum _ {i = 1} ^ n \sum _ {j = 1} ^ m [(i, j) = d]\\
& = \sum _ {d = 1} ^ n d ^ k \sum _ {i = 1} ^ {\lfloor \frac n d \rfloor } \sum _ {j = 1} ^ {\lfloor \frac m d \rfloor } [(i, j) = 1] \\
& = \sum _ {d = 1} ^ n d ^ k \sum _ {i = 1} ^ {\lfloor \frac n d \rfloor } \sum _ {j = 1} ^ {\lfloor \frac m d \rfloor } \sum _ {x | i, x | j} \mu(x)\\
& = \sum _ {d = 1} ^ n d ^ k \sum _ {x = 1} ^ n \mu(x)\lfloor \frac n {xd} \rfloor \lfloor \frac m {xd} \rfloor \\
& = \sum _ {T = 1} ^ n \lfloor \frac n T \rfloor \lfloor \frac m T \rfloor \sum _ {d | T}d ^ k \mu(\frac T d)
\end {aligned}
$$

式子最后是两个积性函数的卷积，依然是积性函数，可以筛出来。

查看代码

#include <cstdio>
#include <algorithm>
using namespace std;
template <class Type>
void read (Type &x)
{
    char c;
    bool flag = false;
    while ((c = getchar()) < '0' || c > '9')
        c == '-' && (flag = true);
    x = c - '0';
    while ((c = getchar()) >= '0' && c <= '9')
        x = (x << 1) + (x << 3) + c - '0';
    flag && (x = ~x + 1);
}
template <class Type, class ...rest>
void read (Type &x, rest &...y) { read(x), read(y...); }
template <class Type>
void write (Type x)
{
    x < 0 && (putchar('-'), x = ~x + 1);
    x > 9 && (write(x / 10), 0);
    putchar('0' + x % 10);
}
typedef long long LL;
const int N = 5e6 + 10, mod = 1e9 + 7;
int n, m, s, cnt, low[N], p[N], f[N];
int binpow (int b, int k = mod - 2)
{
    int res = 1;
    for (; k; k >>= 1, b = (LL)b * b % mod)
        if (k & 1) res = (LL)res * b % mod;
    return res;
}
void init()
{
    low[1] = 1;
    f[1] = 1;
    for (int i = 2; i < N; ++i)
    {
        if (!low[i])
        {
            p[++cnt] = i;
            low[i] = i;
            f[i] = binpow(i, s) - 1;
        }
        for (int j = 1; j <= cnt && i * p[j] < N; ++j)
        {
            if (i % p[j] == 0)
            {
                if (low[i] == i)
                    f[i * p[j]] = f[i] * (f[p[j]] + 1ll) % mod;
                else
                    f[i * p[j]] = (LL)f[i / low[i]] * f[p[j] * low[i]] % mod;
                low[i * p[j]] = low[i] * p[j];
                break;
            }
            f[i * p[j]] = (LL)f[i] * f[p[j]] % mod;
            low[i * p[j]] = p[j];
        }
    }
    for (int i = 2; i < N; ++i)
        (f[i] += f[i - 1]) %= mod;
}
int main ()
{
    int T; read(T, s);
    init();
    while (T--)
    {
        read(n, m);
        if (n > m) swap(n, m);
        int res = 0;
        for (int l = 1, r; l <= n; l = r + 1)
        {
            r = min(n / (n / l), m / (m / l));
            res = (res + (LL)(n / l) * (m / l) % mod * (f[r] - f[l - 1])) % mod;
        }
        write((res + mod) % mod), puts("");
    }
    return 0;
}