BZOJ 2655 calc (组合计数、DP、多项式、拉格朗日插值)

题目链接

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;

}

巴特西

BZOJ 2655 calc (组合计数、DP、多项式、拉格朗日插值)

题目链接

题解

代码

最新文章

热门文章