P6271 [湖北省队互测2014]一个人的数论

P6271 [湖北省队互测2014]一个人的数论

给你 $n=\prod\limits_{i=1}^{w}p_i^{a_i}$，求：

$$\sum_{i=1}^{n}i^d[\gcd(i,n)=1]$$

$1 \leq w \leq 10^3.2 \leq p_i \leq 10^9,1 \leq a_i \leq 10^9$。

先化简原式，有：

$$\begin{aligned} \sum_{i=1}^{n}i^d[\gcd(i,n)=1]&=\sum_{i=1}^{n}i^d\sum_{p|\gcd(i,n)}\mu(p)\\ &=\sum_{i=1}^{n}i^k\sum_{p=1}^{n}\mu(p)[p|i][p|n]\\ &=\sum_{p|n}\mu(p)\sum_{i=1}^{n}i^k[p|i]\\ &=\sum_{p|n}\mu(p)\sum_{i=1}^{\frac{n}{p}}(ip)^d\\ &=\sum_{p|n}\mu(p)k^d\sum_{i=1}^{\frac{n}{p}}i^d \end{aligned}$$

对于 $\sum\limits_{i=1}^{x}i^d$，我们知道它可以表示为一个关于 $x$ 的 $d+1$ 次多项式，设其为：

$$\sum_{i=1}^{d+1}f_ix^i$$

我们可以用高斯消元在 $\mathcal O(d^3)$ 内算出系数 $f_i$ 的值，再算出原式的值。

搞定上面那个后，原式可化为：

$$\begin{aligned} &\sum_{p|n}\mu(p)p^d\sum_{i=0}^{d+1}f_i(\frac{n}{p})^i\\ &=\sum_{i=0}^{d+1}f_i\sum_{p|n}\mu(p)p^d(\frac{n}{p})^i \end{aligned}$$

后面那部分筛一下即可。

#include <bits/stdc++.h>

#define int long long
int read()
{
    int x = 0, f = 1;
    char c = getchar();
    while (c < '0' || c > '9')
    {
        if (c == '-')
            f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9')
    {
        x = x * 10 + c - '0';
        c = getchar();
    }
    return x * f;
}

const int mod = 1e9 + 7;

int qpow(int a, int b)
{
    if (b < 0)
        b += mod - 1;
    int ret = 1;
    for (; b; b >>= 1, a = 1ll * a * a % mod)
        if (b & 1)
            ret = 1ll * ret * a % mod;
    return ret;
}

int inv(int a)
{
    return qpow(a, mod - 2);
}

int m, w;

int p[1007], a[1007];

int mat[105][105], f[105], sum[105];

void init()
{
    for (int i = 1; i <= m + 2; ++i)
    {
        sum[i] = (qpow(i, m) + sum[i - 1]) % mod;
        mat[i - 1][0] = 1, mat[i - 1][m + 2] = sum[i];
        for (int j = 1; j <= m + 1; ++j)
            mat[i - 1][j] = mat[i - 1][j - 1] * i % mod;
    }
}

//行：[0, m + 1]，列：[0, m + 2]

void gauss()
{
    for (int i = 0; i <= m + 1; ++i)
    {
        int r = i;
        for (int j = i + 1; j <= m + 1; ++j)
            if (mat[j][i] > mat[r][i])
                r = j;
        std::swap(mat[r], mat[i]);
        int div = inv(mat[i][i]);
        for (int j = i; j <= m + 2; ++j)
            mat[i][j] = mat[i][j] * div % mod;
        for (int j = i + 1; j <= m + 1; ++j)
        {
            div = mat[j][i];
            for (int k = i; k <= m + 2; ++k)
                mat[j][k] = ((mat[j][k] - mat[i][k] * div % mod) % mod + mod) % mod;
        }
    }
    f[m + 1] = mat[m + 1][m + 2];
    for (int i = m; i >= 0; --i)
    {
        f[i] = mat[i][m + 2];
        for (int j = i + 1; j <= m + 1; ++j)
            f[i] = ((f[i] - mat[i][j] * f[j] % mod) % mod + mod) % mod;
    }
    return;
}

signed main()
{
    m = read(), w = read();
    for (int i = 1; i <= w; ++i)
        p[i] = read(), a[i] = read();
    init();
    gauss();
    int ans = 0;
    for (int i = 0; i <= m + 1; ++i)
    {
        int tmp = 1;
        for (int j = 1; j <= w; ++j)
            tmp = qpow(p[j], a[j] * i % (mod - 1)) * (1ll - qpow(p[j], m - i) + mod) % mod * tmp % mod;
        ans = (ans + tmp * f[i] % mod) % mod;
    }
    printf("%lld\n", ans);
    return 0;
}