2018 Multi-University Training Contest 4 09
题意
给定 \(N\ =\ \Pi_{i\ =\ 1}^t{p_i^{a_i}},\ m\)。令 \(S(i)\ =\ \sum_{j\ =\ 1}^i\ j^m\)。
问 \(\sum_{i\ \leq\ N,\ gcd(i,\ N)\ =\ 1}\ S(i)\) 模 \(998244353\) 的值。
\(1\ \leq\ t\ \leq\ 20,\ 1\ \leq\ m\ \leq\ 10^5\)
做法1
\(\sum_{i\ \leq\ N,\ gcd(i,\ N)\ =\ 1}\ S(i)\ =\ \sum_{d\ |\ N}\ \mu(d)\ \sum_{i\ =\ 1}^{\frac{N}{d}}\ S(id)\)。
令 \(S(x)\ =\ \sum_{j\ =\ 0}^{m\ +\ 1}\ c_j\ x^j,\ B_0\ =\ 1,\ \forall\ i\ >\ 0\ \sum_{i\ =\ 0}^N\ \binom{N\ +\ 1}{i}\ B_i\ =\ 0\)。
\(\sum_{i\ =\ 1}^{N}\ i^k\ =\ \frac{1}{k\ +\ 1}\ \sum_{i\ =\ 0}^k\ B_i\ (N\ +\ 1)^{k\ +\ 1\ -\ i}\ \binom{k\ +\ 1}{i}\)。
\(\sum_{d\ |\ N}\ \mu(d)\ \sum_{i\ =\ 1}^{\frac{N}{d}}\ S(id)\)
\(\ =\ \sum_{d\ |\ N}\ \mu(d)\ \sum_{i\ =\ 1}^{\frac{N}{d}}\ \sum_{j\ =\ 0}^{m\ +\ 1}\ c_j\ (id)^j\)
\(\ =\ \sum_{d\ |\ N}\ \mu(d)\ \sum_{j\ =\ 0}^{m\ +\ 1}\ d^j\ c_j\ \sum_{i\ =\ 1}^{\frac{N}{d}}\ i^j\)
\(\ =\ \sum_{d\ |\ N}\ \mu(d)\ \sum_{j\ =\ 0}^{m\ +\ 1}\ d^j\ c_j\ (\frac{1}{j\ +\ 1}\ \sum_{k\ =\ 0}^j\ B_k\ (\frac{N}{d})^{j\ +\ 1\ -\ k}\ \binom{j\ +\ 1}{k}\ +\ (\frac{N}{d})^j)\)
\(\ =\ \sum_{d\ |\ N}\ \mu(d)\ \sum_{j\ =\ 0}^{m\ +\ 1}\ c_j\ \frac{1}{j\ +\ 1}\ (j\ +\ 1)!\ \sum_{k\ =\ 0}^j\ B_k\ N^{j\ +\ 1\ -\ k}\ d^{k\ -\ 1}\ \frac{1}{k!\ (j\ +\ 1\ -\ k)!}\ +\ \sum_{d\ |\ N}\ \mu(d)\ \sum_{j\ =\ 0}^{m\ +\ 1}\ c_j\ N^j\)
\(\ =\ \sum_{k\ =\ 0}^{m\ +\ 1}\ B_k\ \frac{1}{k!}\ \sum_{j\ =\ k}^{m\ +\ 1}\ c_j\ \frac{1}{j\ +\ 1}\ (j\ +\ 1)!\ N^{j\ +\ 1\ -\ k}\ \frac{1}{(j\ +\ 1\ -\ k)!}\ \sum_{d\ |\ N}\ \mu(d)\ d^{k\ -\ 1}\ +\ \sum_{d\ |\ N}\ \mu(d)\ \sum_{j\ =\ 0}^{m\ +\ 1}\ c_j\ N^j\)
\(\ =\ \sum_{k\ =\ 0}^{m\ +\ 1}\ B_k\ \frac{1}{k!}\ \sum_{j\ =\ k}^{m\ +\ 1}\ c_j\ \frac{1}{j\ +\ 1}\ (j\ +\ 1)!\ N^{j\ +\ 1\ -\ k}\ \frac{1}{(j\ +\ 1\ -\ k)!}\ \Pi_{i\ =\ 1}^t\ (1\ -\ p_i^{k\ -\ 1})\ +\ \sum_{d\ |\ N}\ \mu(d)\ \sum_{j\ =\ 0}^{m\ +\ 1}\ c_j\ N^j\)。
求 \(c_j\) 直接用伯努利数暴力展开 \((x\ +\ 1)^{m\ +\ 1\ -\ i}\) 用fft优化，求伯努利数有 \(B(x)\ =\ \frac{1}{\sum_{i\ \geq\ 0}\ \frac{x^i}{(i\ +\ 1)!}}\)。
代码
#include <bits/stdc++.h>
#ifdef DEBUG
#define debug(...) fprintf(stderr, __VA_ARGS__)
#else
#define debug(...)
#endif
#ifdef __WIN32
#define LLFORMAT "I64"
#else
#define LLFORMAT "ll"
#endif
using namespace std;
const int maxn = (1 << 19) | 10, mod = 998244353, proot = 3, maxt = 22;
inline int pow_mod(int x, int n = mod - 2) { int y = 1; while(n) { if(n & 1) y = (long long) y * x % mod; x = (long long) x * x % mod; n >>= 1; } return y; }
namespace Poly {
struct poly {
vector<int> p;
poly() { p.clear(); }
poly(const vector<int> &q): p(q) {}
poly(int n, int *a) { p.resize(n); for (int i = 0; i < n; ++i) p[i] = a[i]; }
inline int size() const { return p.size(); }
inline void resize(int n) { p.resize(n); return; }
};
void dft(int n, int *a, bool rev) {
for (int i = 0, j = 0; i < n; ++i) {
if(i < j) swap(a[i], a[j]);
for (int k = n >> 1; (j ^= k) < k; k >>= 1);
}
for (int hl = 1, l = 2; l <= n; hl = l, l <<= 1) {
int wn = pow_mod(proot, (mod - 1) / l); if(rev) wn = pow_mod(wn, mod - 2);
for (int i = 0; i < n; i += l) for (int j = 0, *x = a + i, *y = x + hl, w = 1; j < hl; ++j, ++x, ++y, w = (long long) w * wn % mod) {
int t = (long long) *y * w % mod; *y = (*x - t) % mod; *x = (*x + t) % mod;
}
}
if(rev) { int inv = pow_mod(n, mod - 2); for (int i = 0; i < n; ++i) a[i] = (long long) a[i] * inv % mod; }
return;
}
poly operator * (const poly &A, const poly &B) {
static int a[maxn], b[maxn];
int n = A.size(), m = B.size();
if(n < 10 || m < 10 || n + m - 1 < 80) {
int N = n + m - 1; static int c[maxn]; memset(c, 0, sizeof(c[0]) * N);
for (int i = 0; i < n; ++i) for (int x = A.p[i], j = 0; j < m; ++j) c[i + j] = ((long long) x * B.p[j] + c[i + j]) % mod;
return poly(N, c);
}
int N = 1; while(N < n + m - 1) N <<= 1;
for (int i = 0; i < N; ++i) a[i] = i < n ? A.p[i] : 0, b[i] = i < m ? B.p[i] : 0;
dft(N, a, 0); dft(N, b, 0); for (int i = 0; i < N; ++i) a[i] = (long long) a[i] * b[i] % mod; dft(N, a, 1);
return poly(n + m - 1, a);
}
poly operator * (const int &a, const poly &B) {
static int b[maxn]; int n = B.size();
for (int i = 0; i < n; ++i) b[i] = (long long) a * B.p[i] % mod;
return poly(n, b);
}
poly operator - (const poly &A, const poly &B) {
static int a[maxn]; int n = A.size(), m = B.size(), N = max(n, m);
for (int i = 0; i < N; ++i) a[i] = ((i < n ? A.p[i] : 0) - (i < m ? B.p[i] : 0)) % mod;
return poly(N, a);
}
poly operator + (const poly &A, const poly &B) {
static int a[maxn]; int n = A.size(), m = B.size(), N = max(n, m);
for (int i = 0; i < N; ++i) a[i] = ((i < n ? A.p[i] : 0) + (i < m ? B.p[i] : 0)) % mod;
return poly(N, a);
}
poly inv(int n, const poly &A) { // A(x)^-1 mod x^n
if(n == 1) { return poly(vector<int>{pow_mod(A.p[0], mod - 2)}); }
static poly B0, B, tA;
B0 = inv(n + 1 >> 1, A);
if(A.size() < n) tA = A, tA.resize(n);
else tA.p.clear(), tA.p.insert(tA.p.end(), A.p.begin(), A.p.begin() + n);
B = 2 * B0 - B0 * B0 * tA;
B.resize(n);
return B;
}
poly rev(const poly &A) { static poly B; B = A; reverse(B.p.begin(), B.p.end()); return B; }
poly operator / (const poly &A, const poly &B) {
static poly rA, rB, C, D; int n = A.size(), m = B.size();
rA = rev(A); rB = rev(B); C = inv(n - m + 1, rB); D = C * rA; D.resize(n - m + 1);
return rev(D);
}
poly operator % (const poly &A, const poly &B) { static poly D, ret; D = A / B; ret = A - B * D; ret.resize(B.size() - 1); return ret; }
poly Rem(const poly &A, const poly &B, const poly &D) { static poly ret; ret = A - B * D; ret.resize(B.size() - 1); return ret; }
poly dao(const poly &A) {
static int a[maxn]; int n = A.size() - 1;
for (int i = 0; i < n; ++i) a[i] = (long long) A.p[i + 1] * (i + 1) % mod;
return poly(n, a);
}
poly ji(const poly &A) {
static int a[maxn]; int n = A.size();
for (int i = 1; i <= n; ++i) a[i] = (long long) A.p[i - 1] * pow_mod(i, mod - 2) % mod; a[0] = 0;
return poly(n + 1, a);
}
poly ln(int n, const poly &A) { // ln(A(x)) mod x^n
static poly B, C;
C = dao(A); B = inv(n, A);
B = B * C; B = ji(B); B.resize(n);
return B;
}
poly exp(int n, const poly &A) { // e^A(x) mod x^n
if(n == 1) { return poly(vector<int>{1}); }
static poly B0, B, C, tA;
B0 = exp(n + 1 >> 1, A);
if(A.size() < n) tA = A, tA.resize(n);
else tA.p.clear(), tA.p.insert(tA.p.end(), A.p.begin(), A.p.begin() + n);
C = ln(n, B0);
B = B0 * (poly(vector<int>{1}) - C + tA);
B.resize(n);
return B;
}
poly pow_mod(int n, const poly &A, int N) { // A(x)^N mod x^n
static poly B;
int na = A.size(), i = 0;
while(i < na && A.p[i] == 0) ++i;
if(i == na || (long long) N * i >= n) return poly(vector<int>(n, 0));
if(i) {
B.p.clear();
for (int j = i; j < na; ++j) B.p.push_back(A.p[j]);
static poly C;
C = pow_mod(n - N * i, B, N);
B.p = vector<int>(N * i, 0);
B.p.insert(B.p.end(), C.p.begin(), C.p.end());
return B;
}
if(A.p[0] != 1) {
int t = ::pow_mod(A.p[0], mod - 2), s = ::pow_mod(A.p[0], N); B.resize(na);
for (int i = 0; i < na; ++i) B.p[i] = (long long) A.p[i] * t % mod;
B = pow_mod(n, B, N);
for (int i = 0; i < n; ++i) B.p[i] = (long long) s * B.p[i] % mod;
return B;
}
B = ln(n, A);
for (int i = 0; i < n; ++i) B.p[i] = (long long) B.p[i] * N % mod;
return exp(n, B);
}
}
using namespace Poly;
int fac[maxn], ifac[maxn], B[maxn], m, t, a[maxt], p[maxt], c[maxn], h[maxn], N, q[maxt], C[maxt][maxt];
void genfac() {
int n = 100100;
fac[0] = 1;
for (int i = 1; i <= n; ++i) fac[i] = (long long) fac[i - 1] * i % mod;
ifac[n] = pow_mod(fac[n]);
for (int i = n; i; --i) ifac[i - 1] = (long long) ifac[i] * i % mod;
return;
}
void genC() {
int n = maxt - 1;
for (int i = 0; i <= n; ++i) for (int j = C[i][0] = 1; j <= i; ++j) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
return;
}
void genB() {
static poly f;
int n = 100010;
f = poly(n, ifac + 1);
f = inv(n, f);
for (int i = 0; i < n; ++i) B[i] = (long long) f.p[i] * fac[i] % mod;
return;
}
void genc() {
static poly f, g;
int n = m + 2;
f.resize(n); g.resize(n);
for (int i = 0; i < n; ++i) f.p[i] = (long long) B[i] * ifac[i] % mod, g.p[i] = ifac[i];
f = f * g;
for (int i = 0; i <= m + 1; ++i) c[i] = (long long) fac[m] * ifac[i] % mod * f.p[m + 1 - i] % mod;
c[0] = ((long long) c[0] - (long long) fac[m] * ifac[0] % mod * B[m + 1] % mod * ifac[m + 1]) % mod;
return;
}
void genh() {
static poly f, g;
int n = m + 2;
f.resize(n); g.resize(n);
for (int t = N, i = 0; i < n; ++i) f.p[i] = (long long) c[i] * fac[i] % mod, g.p[m + 1 - i] = (long long) t * ifac[i + 1] % mod, t = (long long) t * N % mod;
f = f * g;
for (int i = 0; i <= m + 1; ++i) h[i] = (long long) f.p[i + m + 1] * B[i] % mod * ifac[i] % mod;
return;
}
void solve() {
scanf("%d%d", &m, &t);
N = 1;
for (int i = 1; i <= t; ++i) scanf("%d%d", p + i, a + i), N = (long long) pow_mod(p[i], a[i]) * N % mod, q[i] = pow_mod(p[i]);
genc();
genh();
int ans = 0, sum = 0;
for (int s = 1, i = 0; i <= m + 1; ++i) {
int fuck = h[i];
for (int j = 1; j <= t; ++j) fuck = (long long) fuck * (1 - q[j]) % mod, q[j] = (long long) q[j] * p[j] % mod;
ans = (ans + fuck) % mod;
sum = (sum + (long long) s * c[i]) % mod;
s = (long long) s * N % mod;
}
int smu = 0;
for (int i = 0; i <= t; ++i) {
smu = (((i & 1) ? -1 : 1) * C[t][i] + smu) % mod;
}
ans = (((long long) sum * smu + ans) % mod + mod) % mod;
printf("%d\n", ans);
return;
}
int main() {
int T;
scanf("%d", &T);
genfac();
genC();
genB();
while(T--) solve();
return 0;
}
