题目
给定n（以唯一分解的形式给出）和d，求[1, n]中与n互质的数的d次方和。
结果模1e9+7。
题解
由于n非常大，是以质数幂次方给出。因此应该是要整出一个积性函数，然后拆开求解。
推式子......
[sumlimits_{i=1}^n{i^d[gcd(i, n)=1]}
]
[sumlimits_{i=1}^n{i^dsumlimits_{d_1|gcd(i, j)}{mu(d_1)}}
]
[sumlimits_{d_1|n}{mu(d_1)sumlimits_{i=1}^n{i^d[d_1|i]}}
]
[sumlimits_{d_1|n}{mu(d_1)d_1^dsumlimits_{i=1}^{frac{n}{d_1}}{i^d}}
]
化到这里发现到头了，但(sumlimits_{i=1}^{frac{n}{d_1}}{i^d})不是积性函数。故还要进行进一步的转化。
(sumlimits_{i=1}^{frac{n}{d_1}}{i^d})是d次幂和，可以化为一个关于(frac{n}{d_1})的多项式。即
[sumlimits_{d_1|n}{mu(d_1)d_1^dsumlimits_{i=1}^{k}{a_i(frac{n}{d_1})^i}}
]
[sumlimits_{i=1}^{k}{a_isumlimits_{d_1|n}{mu(d_1)d_1^d(frac{n}{d_1})^i}}
]
现在，(sumlimits_{d_1|n}{mu(d_1)d_1^d(frac{n}{d_1})^i})就是积性函数了，可以愉快地拆开求解。
至于多项式系数(a_i)怎么求，见伯努利数，当然也可以用高斯消元或者拉格朗日插值法。
#include <bits/stdc++.h>
#define endl '
'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define mp make_pair
#define seteps(N) fixed << setprecision(N)
typedef long long ll;
using namespace std;
/*-----------------------------------------------------------------*/
ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f
const int N = 1000 + 10;
const int M = 1e9 + 7;
const double eps = 1e-5;
ll B[N];
ll C[N][N];
ll inv[N];
void init() {
// 预处理组合数
for (int i = 0; i < N; i++) {
C[i][0] = C[i][i] = 1;
for (int k = 1; k < i; k++) {
C[i][k] = (C[i - 1][k] % M + C[i - 1][k - 1] % M) % M;
}
}
// 预处理逆元
inv[1] = 1;
for (int i = 2; i < N; i++) {
inv[i] = (M - M / i) * inv[M % i] % M;
}
// 预处理伯努利数
B[0] = 1;
for (int i = 1; i < N; i++) {
ll ans = 0;
if (i == N - 1) break;
for (int k = 0; k < i; k++) {
ans += C[i + 1][k] * B[k];
ans %= M;
}
ans = (ans * (-inv[i + 1]) % M + M) % M;
B[i] = ans;
}
}
ll aval(int n, int p) {
int k = n + 1 - p;
return inv[n + 1] * C[n + 1][k] % M * B[k] % M;
}
int p[N], a[N];
int rp[N];
inline ll qpow(ll a, ll b, ll m) {
ll res = 1;
while(b) {
if(b & 1) res = (res * a) % m;
a = (a * a) % m;
b = b >> 1;
}
return res;
}
int main() {
IOS;
init();
int d, w;
cin >> d >> w;
ll n = 1;
for(int i = 1; i <= w; i++) {
cin >> p[i] >> a[i];
n = n * qpow(p[i], a[i], M) % M;
rp[i] = qpow(p[i], M - 2, M);
}
ll ans = 0;
ll pn = n;
for(int i = 1; i <= d + 1; i++) {
ll f = 1;
for(int j = 1; j <= w; j++) {
f = f * (1 - qpow(p[j], d, M) % M * qpow(rp[j], i, M) % M) % M;
}
ans += pn * aval(d, i) % M * f % M;
ans %= M;
pn = pn * n % M;
}
cout << (ans + M) % M << endl;
}