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

题意

定义 $f_d(n)=\sum\limits_{i=1}^ni^d\left[\gcd(i,n)=1\right]$，给定 $n=\prod\limits_{i=1}^wp_i^{\alpha_i}$，求 $f_d(n)$，答案对 $10^9+7$ 取模。

数据范围：$1\le d \le 10^2,1\le w\le10^3,1\le p_i,\alpha_i\le10^9,p_i\in\mathbf{Prime}$

题解

$$\begin{aligned}\sum_{i=1}^n i^d\left[\gcd(i,n)=1\right]&=\sum_{i=1}^n i^d\sum_{x|\gcd(i,n)}\mu(x)\\&=\sum_{i=1}^ni^d\sum_{x|i,x|n}\mu(x)\end{aligned}$$ 更换枚举顺序

$$\begin{aligned}\sum_{i=1}^n i^d\left[\gcd(i,n)=1\right]&=\sum\limits_{x|n}\mu(x)\sum\limits_{i=1}^{\left\lfloor{\frac{n}{x}}\right\rfloor}(ix)^d\\&=\sum\limits_{x|n}x^d\mu(x)\sum_{i=1}^{\left\lfloor{\frac{n}{x}}\right\rfloor}i^d\end{aligned}$$

看右半部分这个式子，令 $f(k)=\sum\limits_{i=1}^ki^d$，显然可以用系数表示成一个关于 $k$ 的 $d+1$ 次多项式，设 $f(k)=\sum\limits_{i=0}^{d+1}c_ik^i$，这个多项式的系数可以用高斯消元，lagrange 插值或伯努利数求出，时间复杂度为 $O(d^2)$ 或 $O(d^3)$。

带进去

$$\begin{aligned}\sum_{i=1}^n i^d\left[\gcd(i,n)=1\right]&=\sum\limits_{x|n}x^d\mu(x)f(\left\lfloor{\frac{n}{x}}\right\rfloor)\\&=\sum\limits_{x|n}x^d\mu(x)\sum\limits_{i=0}^{d+1}c_i\left\lfloor{\frac{n}{x}}\right\rfloor^i\\&=\sum\limits_{i=0}^{d+1}c_in^i \sum\limits_{x|n}x^{d-i}\mu(x)\end{aligned}$$

设 $g(n)=\sum\limits_{x|n}x^{d-i}\mu(x)$，这个函数显然是积性的，所以有 $g(n)=\prod\limits_{j=1}^wg(p_j^{\alpha_j})$。

在求 $g(p_j^{\alpha_j})$ 时，根据 $\mu$ 函数的定义，含有完全平方因子时其值为 $0$，所以当且仅当 $x=1$ 或 $x=p_j$ 时，$x^{d-i}\mu(x)$ 对 $g(p_j^{\alpha_j})$ 是有贡献的，即 $g(p_j^{\alpha_j})=1^{d-i}\mu(1)+{p_j}^{d-i}\mu(p_j)=1-p_j^{d-i}$。

所以$$\begin{aligned}\sum_{i=1}^n i^d\left[\gcd(i,n)=1\right]&=\sum\limits_{i=0}^{d+1}c_in^ig(n)\\&=\sum\limits_{i=0}^{d+1}c_in^i\prod\limits_{j=1}^wg(p_j^{\alpha_j})\\&=\sum\limits_{i=0}^{d+1}c_in^i\prod\limits_{j=1}^w(1-p_j^{d-i})\end{aligned}$$

于是可以 $O(d^2+dw\log\text{mod})$ 或 $O(d^3+dw\log\text{mod})$ 求出答案。

#include <bits/stdc++.h>
using namespace std;
namespace IO {
    //read and write
} using namespace IO;
const int N = 1e3 + 10, P = 1e9 + 7;
int qpow(int a, int k) {
    int res = 1;
    while(k) {
        if(k & 1) res = 1ll * res * a % P;
        a = 1ll * a * a % P;
        k >>= 1; 
    }
    return res;
}
int d, w, n, a[N][N], p[N], q[N], c[N], ans;
void build() {
    for(int i = 0; i <= d + 1; i++) {
        for(int j = 0; j <= d + 1; j++) a[i][j] = qpow(i + 1, j);
        for(int j = 1; j <= i + 1; j++)
            a[i][d + 2] = (a[i][d + 2] + qpow(j, d)) % P;
    }
}
void gauss(int n) {
    for(int i = 0, mx, inv; i <= n; i++) {
        mx = i;
        for(int j = i + 1; j <= n; j++)
            if(a[mx][i] < a[j][i]) mx = j;
        swap(a[i], a[mx]);
        inv = qpow(a[i][i], P - 2) % P;
        for(int j = i + 1; j <= n + 1; j++)
            a[i][j] = 1ll * a[i][j] * qpow(a[i][i], P - 2) % P;
        for(int j = 0; j <= n; j++)
            if(i != j)
                for(int k = i + 1; k <= n + 1; k++)
                    a[j][k] = (a[j][k] - 1ll * a[i][k] * a[j][i] % P + P) % P;
    }
    for(int i = 0; i <= n; i++) c[i] = a[i][n + 1]; 
}
int main() {
    d = read(), w = read();
    build(); gauss(d + 1);
    for(int i = 1; i <= w; i++)
        p[i] = read(), q[i] = read(), n = 1ll * n * qpow(p[i], q[i]) % P;
    for(int i = 0, res, t = 1; i <= d + 1; t = 1ll * t * n % P, i++) {
        res = 1ll * c[i] * t % P;
        for(int j = 1; j <= w; j++)
            res = 1ll * res * (1 - qpow(p[j], d - i + P - 1) + P) % P;
        ans = (ans + res) % P;
    }
    printf("%d\n", ans);
    return 0;
}