题目链接

https://www.lydsy.com/JudgeOnline/problem.php?id=2655

题解

据说有一种神仙容斥做法，但我不会。
以及貌似网上大多数人的dp和我的做法都不一样。

下面讲我的做法：
首先由于元素互不相同，那么显然可以先不考虑顺序。
所以要求的就是\(n![x^n]\prod^{m}_{i=1}(1+ix)\) (直接莽上生成函数是不是有点……)
于是发现这个东西和第一类斯特林数生成函数几乎一样，也可以轻易写出递推式\(dp[i][j]=dp[i-1][j]+dp[i-1][j-1]\times i\)

有一个结论是，\(dp[i][j]\)是关于\(i\)的不超过\(2j\)次多项式。
感性理解的话，就是从\(1\)到\(i\)里选\(j\)个，求乘积之和，\(1\)到\(i\)里选\(j\)个一共有\(i\choose j\)种选法，这显然是\(j\)次多项式，再求\(j\)个不超过\(i\)的数的乘积显然也是\(j\)次，那么总共就是\(2j\)次。
于是求出前\(2n\)项，Lagrange插值计算即可。
(所以这其实是一种求第一类斯特林数\(\begin{bmatrix}n\\m\end{bmatrix}\) (\(n-m\)较小)的新方法？)

时间复杂度\(O(n^2)\).

代码

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cassert>
#include<iostream>
#define llong long long
using namespace std;const int N = 1000;
llong fact[N+3],finv[N+3];
llong dp[N+3][N+3];
llong n,m,P;llong quickpow(llong x,llong y)
{llong cur = x,ret = 1ll;for(int i=0; y; i++){if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur%P;}cur = cur*cur%P;}return ret;
}
llong mulinv(llong x) {return quickpow(x,P-2);}namespace Lagrange
{llong ax[N+3],ay[N+3],poly[N+3];llong aux[N+3],aux2[N+3];void lagrange(int n){aux[0] = 1ll;for(int i=0; i<=n; i++){for(int j=i+1; j>0; j--){aux[j] = (aux[j-1]-aux[j]*ax[i]%P+P)%P;}aux[0] = P-aux[0]*ax[i]%P;}for(int i=0; i<=n; i++){llong coe = 1ll;for(int j=0; j<=n; j++){if(i==j) continue;coe = coe*(ax[i]-ax[j]+P)%P;}coe = mulinv(coe);for(int j=0; j<=n+1; j++) aux2[j] = aux[j];for(int j=n; j>=0; j--){poly[j] = (poly[j]+ay[i]*aux2[j+1]%P*coe)%P;aux2[j] = (aux2[j]+aux2[j+1]*ax[i])%P;}}}llong calc(int n,llong x){llong ret = 0ll;for(int i=n; i>=0; i--){ret = (ret*x+poly[i])%P;}return ret;}void clear(int n){for(int i=0; i<=n+1; i++) aux[i] = aux2[i] = poly[i] = 0ll;}
}int main()
{scanf("%lld%lld%lld",&m,&n,&P);fact[0] = 1ll; for(int i=1; i<=N; i++) fact[i] = fact[i-1]*i%P;finv[N] = quickpow(fact[N],P-2); for(int i=N-1; i>=0; i--) finv[i] = finv[i+1]*(i+1)%P;dp[0][0] = 1ll;for(int i=1; i<=n+n; i++){dp[i][0] = 1ll;for(int j=1; j<=i; j++){dp[i][j] = (dp[i-1][j]+dp[i-1][j-1]*i)%P;}}for(int i=0; i<=n+n; i++){Lagrange::ax[i] = i;Lagrange::ay[i] = dp[i][n];}Lagrange::lagrange(n+n);llong ans = Lagrange::calc(n+n,m)*fact[n]%P;printf("%lld\n",ans);return 0;
}