[洛谷P4463] calc （生成函数）

首先注意到题目中 (a) 数组是有序的，那我们只用算有序的方案乘上 (n!) 即可。

而此时的答案显然

[Ans=[x^n](1+x)(1+2x)dots (1+Ax)=prod_{i=1}^A(1+ix) ]

取对数把乘法变加法，即

[prod_{i=1}^A(1+ix)=exp(sum_{i=1}^Aln(1+ix)) ]

这里有 (ln) 的展开式

[-ln(1-x)=sum_{i=1}^inftyfrac{x^i}{i} ]

故有

[ln(1+ix)\=ln(1-(-ix))\=-sum_{k=1}^infty frac{(-ix)^k}{k}\=sum_{k=1}^infty frac{(-1)^{k+1}i^k}{k}x^k ]

则

[sum_{i=1}^A ln(1+ix)\ =sum_{k=1}^infty frac{(-1)^{k+1}sum_{i=1}^Ai^k}{k}x^k ]

自然数幂和可以用某种方法（插值、伯努利数之类）算出来。

最后还要多项式 exp，直接 (O(n^2)) 算。

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N = 505;

int n, M, P, inv[N], Bo[N], C[N][N], a[N], b[N];

typedef vector<int> poly;

poly F[N];

int calc(const poly&a, int x)
{
	int y = 0;
	for(int i = a.size() - 1; i >= 0; i --)
		y = (1LL * y * x + a[i]) % P;
	return y;
}

int main()
{
	scanf("%d%d%d",&M,&n,&P);
	for(int i = 0; i <= n + 1; i ++)
	{
		C[i][0] = 1;
		for(int j = 1; j <= i; j ++)
			C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % P;
	}
	inv[1] = 1;
	for(int i = 2; i <= n + 1; i ++)
		inv[i] = (ll) inv[P % i] * (P - P / i) % P;
	Bo[0] = 1;
	for(int i = 1; i <= n; i ++)
	{
		int t = 0;
		for(int j = 0; j < i; j ++)
			t = (t + (ll)Bo[j] * C[i + 1][j]) % P;
		Bo[i] = (ll)(P - inv[i + 1]) * t % P;
	}
	F[0] = poly{0, 1};
	for(int i = 1; i <= n; i ++)
	{
		F[i].resize(i + 2);
		for(int j = 1; j <= i + 1; j ++)
		{
			F[i][j] = (ll) Bo[i + 1 - j] * C[i + 1][j] % P * inv[i + 1] % P;
			if((i + 1 - j) & 1)
				F[i][j] = (P - F[i][j]) % P;
		}
	}
	for(int i = 1; i <= n; i ++)
	{
		a[i] = (ll)inv[i] * calc(F[i], M) % P;
		if(~i&1) a[i] = (P - a[i]) % P;
	}
	for(int i = 1; i <= n; i ++)
		a[i - 1] = (ll)i * a[i] % P;
	b[0] = 1;
	for(int i = 1; i <= n; i ++)
	{
		for(int j = 0; j < i; j ++)
			b[i] = (b[i] + (ll)b[j] * a[i - j - 1]) % P;
		b[i] = (ll)b[i] * inv[i] % P;
	}
	int ans = b[n];
	for(int i = 2; i <= n; i ++) ans = (ll)ans * i % P;
	printf("%d", ans);
	return 0;
}

注意到复杂度瓶颈在于对所有 (kin [1,n]) 预处理自然数幂和。
这东西的EGF

[sum_{k=0}^infty sum_{i=0}^A i^k frac{x^k}{k!}\ =sum_{i=0}^A sum_{k=0}^infty frac{(ix)^k}{k!}\ = sum_{i=0}^A e^{ix}\= frac{e^{(A+1)x}-1}{e^x-1} ]

多项式求逆即可，整个复杂度也是 (O(nlogn))。
下面是模数 (998244353) 的代码。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<functional>
using namespace std;

typedef long long ll;
const int N = 1 << 20 | 3, P = 998244353;

int fpow(int a, int b)
{
	ll x = 1, o = a;
	for(; b; b >>= 1, o = o * o % P)
		if(b & 1) x = x * o % P;
	return x;
}

int V, n, fac[N], ifac[N], inv[N], a[N], b[N], c[N];

namespace poly
{

int Len, sz, rev[N], w[N];

void prepare(int n)
{
	for(Len = 1, sz = 0; Len <= n; Len <<= 1, sz ++);
	for(int i = 0; i < Len; i ++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (sz - 1));
	int wn = fpow(3, (P - 1) / Len);
	w[Len / 2] = 1;
	for(int i = Len / 2 + 1; i < Len; i ++)
		w[i] = 1LL * w[i - 1] * wn % P;
	for(int i = Len / 2 - 1; i >= 0; i --)
		w[i] = w[i << 1];
}

void DFT(int n, int *a, int T)
{
	static unsigned long long F[N];
	int shift = sz - __builtin_ctz(n), x;
	for(int i = 0; i < n; i ++)
		F[rev[i] >> shift] = a[i];
	for(int l = 1; l < n; l <<= 1)
		for(int i = 0; i < n; i += l << 1)
			for(int j = 0; j < l; j ++)
			{
				x = F[i + j + l] * w[j + l] % P;
				F[i + j + l] = F[i + j] + P - x;
				F[i + j] += x;
			}
	for(int i = 0; i < n; i ++)
		a[i] = F[i] % P;
	if(T)
	{
		x = fpow(n, P - 2);
		for(int i = 0; i < n; i ++)
			a[i] = (ll) a[i] * x % P;
		reverse(a + 1, a + n);
	}
}

void Inverse(int n, int *a, int *b)
{
	if(n == 1)
	{
		*b = fpow(*a, P - 2);
		return;
	}
	Inverse((n + 1) >> 1, a, b);
	static int c[N], len;
	for(len = 1; len < n << 1; len <<= 1);
	for(int i = 0; i < len; i ++)
		i < n ? c[i] = a[i] : c[i] = b[i] = 0;
	DFT(len, b, 0);
	DFT(len, c, 0);
	for(int i = 0; i < len; i ++)
		b[i] = 1LL * b[i] * (P + 2 - 1LL * b[i] * c[i] % P) % P;
	DFT(len, b, 1);
	for(int i = n; i < len; i ++) b[i] = 0;
}

void Exp(int n, int *a, int *b)
{
	static int c[N], d[N];
	for(int i = 1; i < n; i ++)
		c[i - 1] = (ll) a[i] * i % P;
	c[n - 1] = 0;
	for(int i = 0; i < n; i ++)
		d[i] = 0;
	function<void(int,int)> solve = [&](int l, int r)
	{
	    if(l == r)
		{
	        if(!l)
	        	d[l] = 1;
	        else
	        	d[l] = (ll)d[l] * inv[l] % P;
	        return;
	    }
	    int mid = (l + r) / 2;
	    solve(l, mid);
	    static int A[N], B[N];
	    int L = 1;
	    while(L <= r - l) L <<= 1;
	    memset(A, 0, L << 2);
	    memset(B, 0, L << 2);
	    memcpy(A, d + l, (mid - l + 1) << 2);
	    memcpy(B, c, (r - l) << 2);
	    DFT(L, A, 0);
	    DFT(L, B, 0);
	    for(int i = 0; i < L; i ++)
	    	A[i] = (ll) A[i] * B[i] % P;
	    DFT(L, A, 1);
	    for(int i = mid + 1; i <= r; i ++)
	    	d[i] = (d[i] + A[i - l - 1]) % P;
	    solve(mid + 1, r);
	};
	solve(0, n - 1);
	memcpy(b, d, n << 2);
}

}

int main()
{
	scanf("%d %d", &V, &n);
	poly::prepare(n + n);
	fac[0] = ifac[0] = 1;
	for(int i = 1; i <= n + 1; i ++)
	{
		fac[i] = (ll)fac[i - 1] * i % P;
		inv[i] = (i != 1 ?
			(ll)inv[P % i] * (P - P / i) % P : 1);
		ifac[i] = (ll)ifac[i - 1] * inv[i] % P;
	}
	// (e^x)^(V+1)-1 / e^x-1
	for(int i = 0, e = V + 1; i <= n; i ++, e = 1LL * e * (V + 1) % P)
	{
		a[i] = 1LL * e * ifac[i + 1] % P;
		b[i] = ifac[i + 1];
	}
	poly::Inverse(n + 1, b, c);
	int L = poly::Len;
	poly::DFT(L, a, 0);
	poly::DFT(L, c, 0);
	for(int i = 0; i != L; i ++) a[i] = 1LL * a[i] * c[i] % P;
	poly::DFT(L, a, 1);
	a[0] = 0;
	for(int i = 1; i <= n; i ++)
	{
		a[i] = (ll)a[i] * fac[i] % P * inv[i] % P;
		if((i + 1) & 1) a[i] = (P - a[i]) % P;
	}
	poly::Exp(n + 1, a, b);
	for(int i = 1; i <= n; i ++)
		printf("%d
", int(1LL * b[i] * fac[i] % P));
	return 0;
}