先了解一个柿子：(inom{n}{k} imes k^{underline m} = inom{n-m}{k-m} imes n^{underline m})。
证明：
[egin{aligned}
inom{n}{k} imes k^{underline m} &= frac{n!}{k!(n-k)!} imes frac{k!}{(k-m)!} \
&= frac{n!}{(n-k)!(k-m)!} \
&= frac{n!(n-m)!}{(n-k)!(k-m)!(n-m)!} \
&= inom{n-m}{k-m}frac{n!}{(n-m)!} \
&= inom{n-m}{k-m}n^{underline m}
end{aligned}]
我们发现题目中的多项式不是很好算，于是把它改成下降幂形式：
[f(k)=sumlimits_{i=1}^m a_ik^i=sumlimits_{i=1}^m b_ik^{underline i}
]
考虑怎么求 (b_i)：
[egin{aligned}
sumlimits_{i=0}^m a_ik^i &= sumlimits_{i=0}^m a_isumlimits_{j=0}^i inom{k}{j}{irace j}j! \
&= sumlimits_{i=0}^m a_isumlimits_{j=0}^i k^{underline j}{irace j}\
&= sumlimits_{j=0}^m k^{underline j}sumlimits_{i=j}^m {irace j}a_i
end{aligned}]
所以 (b_i=sumlimits_{j=i}^m {jrace i}a_j)。
然后就开始推式子：
[egin{aligned}
sumlimits_{k=0}^nsumlimits_{i=0}^m b_ik^{underline{i}}x^kinom{n}{k} &= sumlimits_{k=0}^nsumlimits_{i=0}^m inom{n-i}{k-i}n^{underline{i}}b_ix^k \
&= sumlimits_{i=0}^m b_in^{underline{i}}sumlimits_{k=0}^n inom{n-i}{k-i}x^k\
&= sumlimits_{i=0}^m b_in^{underline{i}}sumlimits_{k=0}^{n-i}x^{k+i}inom{n-i}{i} \
&= sumlimits_{i=0}^m b_in^{underline{i}} x^isumlimits_{k=0}^{n-i}inom{n-i}{i}x^k1^{n-i-k} \
&= sumlimits_{i=0}^m b_in^{underline{i}}x^i(x+1)^{n-i}
end{aligned}]
可以 (mathcal{O}(m^2)) 求解。
代码：
#include <bits/stdc++.h>
#define DEBUG fprintf(stderr, "Passing [%s] line %d
", __FUNCTION__, __LINE__)
#define File(x) freopen(x".in","r",stdin); freopen(x".out","w",stdout)
using namespace std;
typedef long long LL;
typedef pair <int, int> PII;
typedef pair <int, PII> PIII;
template <typename T>
inline T gi()
{
T f = 1, x = 0; char c = getchar();
while (c < '0' || c > '9') {if (c == '-') f = -1; c = getchar();}
while (c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
return f * x;
}
const int INF = 0x3f3f3f3f, N = 1003, M = N << 1;
int n, x, p, m;
LL a[N], b[N], s[N][N], ans;
inline LL qpow(LL g, LL u)
{
LL res = 1;
while (u)
{
if (u & 1) res = res * g % p;
g = g * g % p, u >>= 1;
}
return res;
}
int main()
{
//File("");
n = gi <int> (), x = gi <int> (), p = gi <int> (), m = gi <int> ();
s[0][0] = 1;
for (int i = 1; i <= m; i+=1)
for (int j = 1; j <= i; j+=1)
s[i][j] = (s[i - 1][j - 1] + 1ll * j * s[i - 1][j] % p) % p;
for (int i = 0; i <= m; i+=1) a[i] = gi <int> ();
for (int i = 0; i <= m; i+=1)
for (int j = i; j <= m; j+=1)
b[i] = (b[i] + s[j][i] * a[j] % p) % p;
for (int i = 0; i <= m; i+=1)
{
LL tmp = 1;
for (int j = n - i + 1; j <= n; j+=1) tmp = tmp * j % p;
LL now = b[i] % p * tmp % p * qpow(x, i) % p * qpow(x + 1, n - i) % p;
ans = (ans + now) % p;
}
printf("%lld
", ans);
return 0;
}