[bzoj3601]一个人的数论【高斯消元】【莫比乌斯反演】

【题目链接】
　　https://www.lydsy.com/JudgeOnline/problem.php?id=3601
【题解】
　　首先给定的式子不是一个积性函数（一开始想都没想写了一发，直接gg）。
　　我们要求的是:$\sum_{i=1}^{n}[gcd(i,n)==1]i^d$
　　看到gcd先反演一下，变成$\sum_{i|n}\mu(i)i^d\sum_{j=1}^{n/i}j^d$
　　考虑后面的和式，一定是一个$d+1$次多项式，可以把这个多项式通过高斯消元求出。
　　记多项式的系数为$A_{0}..A_{d+1}$
　　于是式子变为：$\sum_{i|n}\mu(i)i^d\sum_{j=0}^{d+1}A_{j}(n/i)^j$
　　把$A$提出：$\sum_{j=0}^{d+1}A_{j}\sum_{i|n}\mu(i)i^d(n/i)^j$
　　显然$\mu(i)i^d$是积性函数，$i^j$是积性函数。
　　所以后面的和式可以看做积性函数的卷积，也是积性函数，所以每一个$p$可以分开求。
　　由于$\mu(p^i)==0(i>1)$，所以只要枚举两项。
　　时间复杂度$O(D^3+N*D)$
【代码】

/* - - - - - - - - - - - - - - -
    User :      VanishD
    problem :   [bzoj3601]
    Points :    gauss
- - - - - - - - - - - - - - - */
# include <bits/stdc++.h>
# define    ll      long long
# define    N       110
using namespace std;
const int P = 1e9 + 7;
const int inf = 0x3f3f3f3f, INF = 0x7fffffff;
const ll  infll = 0x3f3f3f3f3f3f3f3fll, INFll = 0x7fffffffffffffffll;
int x[N], y[N], d, w, num, sum[N], A[N][N], ans;
int read(){
    int tmp = 0, fh = 1; char ch = getchar();
    while (ch < '0' || ch > '9'){ if (ch == '-') fh = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9'){ tmp = tmp * 10 + ch - '0'; ch = getchar(); }
    return tmp * fh;
}
int power(int x, int y){
    int i = x; x = 1;
    while (y > 0){
        if (y % 2 == 1) x = 1ll * x * i % P;
        i = 1ll * i * i % P;
        y /= 2;
    }
    return x;
}
void gauss(int n){
    for (int i = 1; i <= n; i++){
        int now, tmp;
        for (int j = i; j <= n; j++)
            if (A[j][i] != 0) now = j;
        if (now != i) swap(A[now], A[i]);
        for (int j = 1; j <= n; j++){
            if (A[j][i] == 0 || j == i) continue;
            tmp = 1ll * A[j][i] * power(A[i][i], P - 2) % P;
            for (int k = i; k <= n + 1; k++)
                A[j][k] = (A[j][k] - 1ll * A[i][k] * tmp % P + P) % P;
        }
    }
    for (int i = 1; i <= n; i++) sum[i - 1] = 1ll * A[i][n + 1] * power(A[i][i], P - 2) % P;
}
int main(){
//  freopen("bzoj3601.in", "r", stdin);
//  freopen("bzoj3601.out", "w", stdout);
    d = read(), w = read();
    num = d + 2;
    for (int i = 1; i <= num; i++){
        x[i] = i, y[i] = (y[i - 1] + power(i, d)) % P;
        A[i][num + 1] = y[i];
        for (int j = 1; j <= num; j++) A[i][j] = power(x[i], j - 1);
    }
    gauss(num);
    for (int i = 1; i <= w; i++){
        int p = read(), a = read();
        for (int j = 0; j <= d + 1; j++)
            sum[j] = 1ll * sum[j] * (1 * power(power(p, a), j) + -1ll * power(p, d) * power(power(p, a - 1), j) % P) % P;
    }
    for (int i = 0; i <= d + 1; i++) ans = (ans + sum[i]) % P;
    printf("%d\n", (ans + P) % P);
    return 0;
}