题目大意：对于每个$k\in[1,t]$，求：
$$
\dfrac{\sum\limits_{i=1}^n\sum\limits_{j=1}^m(a_i+b_j)^k}{nm}
$$
$n,m,t\leqslant10^5$
题解：发现这个$nm$是一个定值，可以先不考虑，先对每一个$k$来求
$$
\begin{align*}
&\sum\limits_{i=1}^n\sum\limits_{j=1}^m(a_i+b_j)^k\\
=&\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{s=0}^k\binom ksa_i^sb_j^{k-s}\\
=&\sum\limits_{s=0}^k\binom ks\sum\limits_{i=1}^na_i^s\sum\limits_{j=1}^mb_j^{k-s}\\
=&\sum\limits_{s=0}^k\dfrac{k!}{s!(k-s!)}\sum\limits_{i=1}^na_i^s\sum\limits_{j=1}^mb_j^{k-s}\\
=&k!\sum\limits_{s=0}^k\sum\limits_{i=1}^n\dfrac{a_i^s}{s!}\sum\limits_{j=1}^m\dfrac{b_j^{k-s}}{(k-s)!}\\
\end{align*}
$$
令$A_k=\sum\limits_{i=0}^n\dfrac{a_i^k}{k!},B_k=\sum\limits_{i=0}^m\dfrac{b_i^k}{k!}$式子就变成了
$$
k!\sum\limits_{s=0}^kA_sB_{k-s}\\
$$
就可以很简单的求出来了，现在的问题变成了如何求$A,B$。接下来就讲求$A$的方法，$B$的相同。
令$A(x)$为$A$的生成函数，即$A(x)=\sum\limits_{k=0}^{\infty}\dfrac{\sum\limits_{i=0}^na_i^k}{k!}x^k$，这一个$k!$也很好解决，下面就省略不写
$$
\begin{align*}
f(x)=&\sum\limits_{k=0}^{\infty}\sum\limits_{i=0}^na_i^kx^k\\
=&\sum\limits_{i=0}^n\sum\limits_{k=0}^{\infty}a_i^kx^k\\
=&\sum\limits_{i=0}^n\dfrac1{1-a_i}
\end{align*}
$$
发现$\ln‘(x)=\dfrac1x,\ln‘(1-ax)=\dfrac{-a_i}{1-a_i}$
令$g(x)=\sum\limits_{i=0}^n\dfrac{-a_i}{1-a_i}$，然后$f(x)=-g(x)x+n$，而$g(x)$可以比较简单的求出来：
$$
\begin{align*}
g(x)=&\sum\limits_{i=0}^n\dfrac{-a_i}{1-a_i}\\
=&\sum\limits_{i=0}^n\ln‘(1-a_ix)\\
=&\ln‘(\prod\limits_{i=1}^n(1-a_ix))
\end{align*}
$$
令$M(x)=\prod\limits_{i=1}^n(1-a_ix)$，这可以用分治$FFT$求出
$$
\begin{align*}
g(x)=&\ln‘(M(x))\\
=&(\int\dfrac{M‘(x)}{M(x)})‘\\
=&\dfrac{M‘(x)}{M(X)}
\end{align*}
$$
然后把$g(x)$转成$f(x)$，再变成指数型生成函数卷一下就好了，注意别忘记上面省略掉的一些东西
卡点：求导写挂，一处地方不开$C++11$毁青春
C++ Code：
#include <algorithm>
#include <cstdio>
#include <iostream>
#include <vector>
#define maxn 262144
const int mod = 998244353;
#define mul(x, y) static_cast<long long> (x) * (y) % mod
namespace Math {
inline int pw(int base, int p) {
static int res;
for (res = 1; p; p >>= 1, base = mul(base, base)) if (p & 1) res = mul(res, base);
return res;
}
inline int inv(int x) { return pw(x, mod - 2); }
}
inline void reduce(int &x) { x += x >> 31 & mod; }
inline void clear(register int *l, const int *r) {
if (l >= r) return ;
while (l != r) *l++ = 0;
}
int inv[maxn], ifac[maxn], fac[maxn];
int A[maxn], B[maxn];
namespace Poly {
#define N maxn
int lim, s, rev[N], Wn[N];
inline void init(const int n) {
lim = 1, s = -1; while (lim < n) lim <<= 1, ++s;
for (register int i = 1; i < lim; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << s;
const int t = Math::pw(3, (mod - 1) / lim);
*Wn = 1; for (register int *i = Wn + 1; i != Wn + lim; ++i) *i = mul(*(i - 1), t);
}
inline void FFT(int *A, const int op = 1) {
for (register int i = 1; i < lim; ++i) if (i < rev[i]) std::swap(A[i], A[rev[i]]);
for (register int mid = 1; mid < lim; mid <<= 1) {
const int t = lim / mid >> 1;
for (register int i = 0; i < lim; i += mid << 1)
for (register int j = 0; j < mid; ++j) {
const int X = A[i + j], Y = mul(A[i + j + mid], Wn[t * j]);
reduce(A[i + j] += Y - mod), reduce(A[i + j + mid] = X - Y);
}
}
if (!op) {
const int ilim = Math::inv(lim);
for (register int *i = A; i != A + lim; ++i) *i = mul(*i, ilim);
std::reverse(A + 1, A + lim);
}
}
std::vector<int> P1[maxn], P2[maxn];
#define __DC_FFT(x, A) void DC_FFT##x(int rt, int l, int r) { if (l == r) { P##x[rt] = {1, A[l]}; return ; } static int C[N], D[N]; const int mid = l + r >> 1, L = rt << 1, R = rt << 1 | 1; DC_FFT##x(L, l, mid), DC_FFT##x(R, mid + 1, r); const int n = P##x[L].size(), m = P##x[R].size(); init(n + m); std::copy(P##x[L].begin(), P##x[L].end(), C), clear(C + n, C + lim); std::copy(P##x[R].begin(), P##x[R].end(), D), clear(D + m, D + lim); FFT(C), FFT(D); for (int i = 0; i < lim; ++i) C[i] = mul(C[i], D[i]); FFT(C, 0); P##x[rt].assign(C, C + n + m - 1); }
__DC_FFT(1, A)
__DC_FFT(2, B)
inline void DER(int *A, int *B, int n) {
B[n - 1] = 0; for (int i = 1; i < n; ++i) B[i - 1] = mul(A[i], i);
}
void INV(int *A, int *B, int n) {
if (n == 1) { *B = Math::inv(*A); return ; }
static int C[N], D[N];
const int len = n + 1 >> 1;
INV(A, B, len), init(len * 3);
std::copy(A, A + n, C), clear(C + n, C + lim);
std::copy(B, B + len, D), clear(D + len, D + lim);
FFT(C), FFT(D);
for (int i = 0; i < lim; ++i) D[i] = (2 - mul(C[i], D[i]) + mod) * D[i] % mod;
FFT(D, 0); std::copy(D + len, D + n, B + len);
}
void solve(std::vector<int> P, int *A, int t) {
static int B[N], C[N];
const int n = P.size();
std::copy(P.begin(), P.end(), A); clear(A + n, A + t);
DER(A, B, n), INV(A, C, t);
init(n + t);
clear(B + n, B + lim), clear(C + t, C + lim);
FFT(B), FFT(C);
for (int i = 0; i < lim; ++i) A[i] = mul(B[i], C[i]);
FFT(A, 0);
for (int i = t; i; --i) reduce(A[i] = -A[i - 1]);
A[0] = n - 1;
for (int i = 0; i < t; ++i) A[i] = mul(A[i], ifac[i]);
}
void TANG_Yx(int n, int m, int t) {
DC_FFT1(1, 0, n - 1), DC_FFT2(1, 0, m - 1);
static int A[N], B[N];
solve(P1[1], A, t), solve(P2[1], B, t);
init(t << 1);
clear(A + t, A + lim), clear(B + t, B + lim);
FFT(A), FFT(B);
for (int i = 0; i < lim; ++i) A[i] = mul(A[i], B[i]);
FFT(A, 0);
const int nm = mul(Math::inv(n), Math::inv(m));
for (int i = 1; i < t; ++i) printf("%lld\n", mul(A[i], fac[i]) * nm % mod);
}
#undef N
}
int n, m, t;
int main() {
std::ios::sync_with_stdio(false), std::cin.tie(0), std::cout.tie(0);
inv[0] = inv[1] = ifac[0] = ifac[1] = fac[0] = fac[1] = 1;
for (int i = 2; i < maxn; ++i) {
inv[i] = mul(mod - mod / i, inv[mod % i]);
ifac[i] = mul(ifac[i - 1], inv[i]);
fac[i] = mul(fac[i - 1], i);
}
std::cin >> n >> m;
for (int i = 0; i < n; ++i) std::cin >> A[i], reduce(A[i] = -A[i]);
for (int i = 0; i < m; ++i) std::cin >> B[i], reduce(B[i] = -B[i]);
std::cin >> t; ++t;
Poly::TANG_Yx(n, m, t);
return 0;
}
原文：https://www.cnblogs.com/Memory-of-winter/p/10372741.html