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;
}

转载于:https://www.cnblogs.com/suncongbo/p/11291510.html